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Abstract 


In recent years, fiber-reinforced materials have risen to prominence in a variety of ap- 
plications spanning from aircraft and automobile engineering to medical technology and 
biological materials. Despite the scientific progress in understanding this class of mate- 
rials achieved in diverse research contributions, the precise modeling and simulation of 
macroscopic systems with dedicated microstructure up to ultimate fracture remains an 
open research field in mechanics. 

The present work presents a novel numerical approach to analyze the thermomechanical 
behavior within composite materials including the inelastic regime up to final failure. 
Therefore, a second-gradient theory is combined with phase-field methods to fracture. 
In particular, we assume that the polymeric matrix material undergoes ductile fracture, 
whereas continuously embedded fibers undergo brittle fracture as it is typical e.g. for 
roving glass reinforced thermoplastics. A hybrid phase-field approach is developed and 
applied in conjunction with a modified Gurson-Tvergaard-Needleman GTN-type plastic- 
ity model accounting for a temperature-dependent growth of voids on microscale. The 
mechanical response of the arising microstructure of the woven fabric gives rise to addi- 
tional higher-order terms, representing homogenized bending contributions of the fibers. 
Eventually, a series of tests is conducted for this physically comprehensive multifield for- 
mulation to investigate various types and sequences of failure within long fiber-reinforced 
polymers. 


Kurzfassung 


In den letzten Jahren haben faserverstárkte Werkstoffe in einer Vielzahl von Anwendun- 
gen vom Flugzeug- und Automobilbau tiber die Medizintechnik bis hin zu biologischen 
Materialien eine immense Bedeutung gewonnen. Trotz der wissenschaftlichen Fortschritte 
im Verstándnis dieser Materialklassen, die in diversen Forschungsbeitrágen erzielt wur- 
den, bleibt die prázise Modellierung und Simulation makroskopischer Systeme mit dedi- 
zierter Mikrostruktur bis hin zum endgiiltigen Bruch ein offenes Forschungsfeld in der 
Mechanik. Die vorliegende Arbeit stellt einen neuartigen numerischen Ansatz vor, um 
das thermomechanische Verhalten in Verbundwerkstoffen einschließlich des inelastischen 
Regimes bis zum endgúltigen Versagen zu analysieren. Dazu wird eine zweite Gradi- 
enten Theorie mit Phasenfeldmethoden fiir den Bruch kombiniert. Insbesondere wird 
davon ausgegangen, dass das polymere Matrixmaterial einen duktilen Bruch erfáhrt, 
während kontinuierlich eingebettete Fasern einen Sprödbruch erfahren, wie er z.B. fiir 
glasfaserverstärkte Thermoplaste typisch ist. Es wird ein hybrider Phasenfeldansatz en- 
twickelt und in Verbindung mit einem modifizierten Gurson-Tvergaard-Needleman GTN- 
Plastizitätsmodell angewendet, das ein temperaturabhängiges Wachstum von Hohlräu- 
men auf der Mikroskala berücksichtigt. Die mechanische Charakterisierung der entste- 
henden Mikrostruktur des Gewebes führt zu zusätzlichen Termen höherer Ordnung, die 
homogenisierte Biegebeiträge der Fasern darstellen. Schließlich wird für diese physikalisch 
umfassende Mehrfeldformulierung eine Versuchsreihe durchgeführt, um verschiedene Ver- 
sagensarten und -folgen in langfaserverstärkten Polymeren zu untersuchen. 
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1 Introduction 


Composite materials have played a significant part in human development throughout 
history. Thus, the combination of clay and straw for the production of bricks is known 
from the excavation of the city of Çatalhöyük [77], the first documented city, dated back 
to about 8000 BC. Whilst straw improves the tensile capacity of clay, it will also act to 
assist in keeping the moisture content of the material more constant throughout, leading 
to less drying shrinkage cracks. 


In the development of polymeric fiber composites, the low density and high strength fiber 
material, typically glass or carbon, is combined with a polymer matrix in order to create 
geometrically defined components. The matrix serves not only to shape the material, 
but also to transfer loads to the fibers and to protect the fibers from environmental 
influences and abrasive wear caused by mutual contact. An early example of the use of 
polymer matrix composites are the Micarta aircraft propellers from around 1920 onwards 
[71]. Due to their excellent cost-performance ratio as well as their adjustable mechanical 
properties, fiber-reinforced composites are nowadays used in a wide range of industrial 
applications. In recent years, new types of additive manufacturing machines have been 
developed and brought to market that are capable of designing customized composites 
by introducing fibers into the produced parts during the manufacturing process. Until 
now, however, the design of structural composite components has been governed by an 
incomplete knowledge, which is mainly based on empiricism. 


However, the dramatic increase in computing power over the past decades has made possi- 
ble more and more realistic simulations of complex materials. We can now resolve sophis- 
ticated micro structures in a multi-physical environment using brute-force computations 
within a continuum mechanical framework. New challenges arise for the implementation 
of composite materials in a suitable simulation environment in order to obtain realistic 
predictions of the material behavior. In particular, the distribution and orientation of 
the fibers on the microscale have a significant influence on the macroscale mechanics. 
Although of great interest to the industry as a whole, numerical prediction of thermome- 
chanical macroscopic systems with dedicated microstructure behavior including arbitrary 
three-dimensional fracture patterns of fiber-reinforced polymers is currently subject to 
large uncertainties and still under research. This is a dramatic information deficit, since 
fiber-reinforced polymers are already in use in vital structures. 


2 1 Introduction 


1.1 State of the art and objectives 


1.1.1 Strain-gradient formulations for fiber-reinforced materials 


Fiber-reinforced materials can be understood as having detailed microstructure because 
the length scales of the embedded fibers are small compared to the surrounding continuum. 
A very general framework for including the effects of microstructures within a continuum 
formulation was presented in the seminal and fundamental work on higher-order theories 
by Mindlin |89, 90], see also the work of Germain [53], Toupin |121, 122], and Eringen 
[21]. 


A number of publications establish a connection between the microstructure and the 
macroscopic formulation as generalized continua, see, among others |81, 49, 60, 105], 
see also |72, 67, 97, 98, 129] for more details on the mathematical structure. Specific 
mechanical problems such as elastic meshes have been treated by Steigmann et al. |118, 
115, 117, 116], for application to panthographic structures see dell'Isola et al. |31, 32]. 


Within the present work a model of a woven fabric as presented in Steigmann |116] will 
be first be embedded into the Kirchhoff-Love shell theory and then be adapted to three- 
dimensional continua. Therein, the in-plane flexural resistance of the fibers is taken into 
account in addition to first-gradient anisotropic effects. We will demonstrate on the basis 
of experimental results that a classical Cauchy continuum theory without higher-gradient 
contributions can only be adapted to a specific fiber orientation and load case, but never 
independently on the orientation. In contrast, the proposed formulation as generalized 
continuum with higher-gradient contributions allows for an independent modeling without 
recalibration of the material for the specific fiber direction. Further discussions on higher- 
order contributions for the constitutive modeling of composites have been presented in 
Asmanoglo and Menzel |10, 11], using the framework as provided in the preliminary work 
of Spencer and Soldatos |114| and Soldatos |113]. Note that we will show here, that 
the formulation as proposed by dell'Isola et al. [32] can be recast into the formulation of 
Asmanoglo and Menzel |10, 11]. 


To obtain the corresponding constitutive information, classical homogenization methods 
are often applied, which are usually limited to the first gradient anisotropic contributions. 
In particular, full-field homogenizations are usually based on a unit cell discretized with 
finite elements (see, e.g., in Fish |46] for a comprehensive review of multiscale methods) 
and have been applied to composites, including in Wang et al. [125], see also the references 
in Mishnaevsky |93]. Basic ideas on the homogenization of continua with micromorphic 
mesostructure can be found in Hirschberger |59]. Recent studies on the homogenization 
of continua with micromorphic properties can be found in Hutter [65]. 


For the spatial discretization, B-spline based shape functions are used, which originate 
from the field of Computer Aided Design (CAD). The corresponding finite element frame- 
work is commonly known as isogeometric analysis (IGA), see Hughes et al. [64] and Cot- 
trell et al. [26]. IGA allows to control the differentiability of the shape functions and is 
therefore suitable for higher-gradient continua, see Fischer et al. |45] for applications to 
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strain-gradient elasticity. Kirchhoff-Love shell formulations have also been analyzed in 
recent years using the concept of IGA, see Kiendl et al. |69, 70], among others. 


1.1.2 Phase-field approach to fracture 


Along this line, Choo and Sun |22] have developed a framework that combines a phase- 
field approach to fracture with a pressure-sensitive plasticity for modeling brittle fracture 
to ductile flow in geological materials. A coupled crystal-plasticity-phase-field fracture 
formulation based on experimental observations was presented in Diehl et al. [33]. A 
gradient-extended plasticity model for overcoming mesh sensitivities, such as that pre- 
sented in De Borst and Mühlhaus [30], Fleck et al. [47], Menzel and Steinmann |82], Geers 
[52], Anand et al. [9], Dimitrijevic and Hackl |34], Miche et al. [87], Aldakheel |3, 5], Wulf- 
inghoff and Böhlke |127], and Forest |48] was applied to phase-field breaking problems in 
Aldakheel et al. [4], Miehe et al. |84], Aldakheel [2], and Dittmann et al. [35]. 


A variety of phenomenological and micromechanical approaches exist in the literature for 
modeling the diverse damage mechanisms in the polymers used here as matrix materials. 
To describe such phenomena, material behavior at the microscale must be included in 
the continuum formulation. In particular, the growth of microvoids prior to ultimate 
fracture must be considered at the macroscale, as rooted in the seminal work of Gurson 
|56, 57]. In this work, a macroscopic flow surface was developed by homogenizing a 
porous representative volume element with assumed rigid plastic flow, which degrades 
with increasing void fraction. Later, this model was used by Tvergaard |123], Tvergaard 
and Needleman [124], Needleman and Tvergaard |96], Leblond et al. [74], Nahshon and 
Hutchinson [95], Xue et al. [128], Li et al. [76], and Huespe et al. [62] to account for damage 
growth, extending the flow criterion function by introducing new material parameters to 
account for nucleation and coalescence effects. In Hutter et al. [63] and Reusch et al. 
|103, 104], nonlocal Gurson models were presented to overcome the nonphysical mesh 
sensitivity in the softening materials. Recent works on the application of this model 
to polyamide, as commonly used for fiber-reinforced thermoplastics, can be found, for 
example, in Selles et al. [109] and Cayzac et al. |20]. Thermomechanical extensions of 
the Gurson-type model with brittle crack propagation in thermoelastic solids were shown 
in Dittmann et al. [39] and Miehe et al. [86]. An extension towards finite-strain thermo- 
poros-plasticity based on the phase-field approach was recently developed in Dittmann et 


al. [36]. 


In this contribution, we introduce a novel hybrid phase-field model to simulate failure 
mechanisms in fiber-reinforced materials. In particular, we assume that the polymeric ma- 
trix material undergoes ductile fracture, whereas continuously embedded fibers undergo 
brittle fracture. The matrix and fiber material are associated with indepedent phase- 
fields. A modified Gurson-Tvergaard-Needleman GTN-type plasticity model accounting 
for a temperature-dependent growth of voids on microscale is furthermore applied to the 
polymeric matrix material. Eventually, a series of tests is conducted for this physically 
comprehensive multifield formulation to investigate different kinds and sequences of failure 
within long fiber-reinforced polymers. 


4 1 Introduction 


1.2 Outline of the thesis 


The outline of the thesis is as follows. Chapter 2 introduces the basic notations and 
concepts of classical continuum mechanics used in this thesis. The approach to sim- 
ulate porous-ductile fracture in isotropic thermo-elasto-plastic solids subjected to large 
deformations is then presented in Chapter 3. Chapter 4 provides an introduction to gen- 
eralized continuum mechanics for microstructured materials. A very specific generalized 
continuum formulation for fiber materials is then outlined in Chapter 5. The complete 
framework for modeling and simulation of fiber-reinforced polymers, including the inelas- 
tic regime up to ultimate failure, is presented in Chapter 6. Finally, conclusions are drawn 
in Chapter 7. 


2 Fundamentals 


This chapter introduces the basic notations and concepts of this work. After recalling 
some tensor operations, the principle equations of non-linear continuum mechanics are 
presented. Next, the kinematics of Kirchhoff-Love shells are outlined. Finally, the phase- 
field approach to non-linear brittle fracture and the isogeometric concept for finite element 
analysis are introduced. 


2.1 Tensor notations 


In order to provide a clear representation of mathematical operations used therein we 
recall and define some specific tensor operations. Let a € R? be a vector, A € R? x R? a 
second-order tensor and 21 € R? x R? x R? a third-order tensor. 


We denote the symmetric, skew-symmetric, spheric, and deviatoric part of a second-order 
tensor, respectively, as 


where I is the identity matrix. 


In addition, we define the gradient with respect to the reference and current configuration 
V(e) and V..(e) of a vector field as 


_ dla]; _ Ola); 
[Va]; = IX], and [Va]; = I]; (2.1) 
and of a second-order tensor field as 
_ Air OLA]; 
[VA] = IK and [Ve Alix = Alas (2.2) 


respectively. The divergence operator with respect to the reference configuration Div(e) 
or V - (e) of a second-order tensor field A and third-order tensor field A is defined as 


OLA]: OW li. 


ee AX 


(2.3) 


and [V 5 Al. = 
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respectively. Finally, the double contractions of a second-order and third-order tensor, 
i.e. a = A: Y and b= 2%: A, are defined as 


la], = [A]; : 2Qije and [bl = [Wize : LA] yx, (2.4) 


respectively. 


2.2 Continuum mechanics 


This section recalls the fundamental concepts of classical continuum mechanics at finite 
strains, also referred to as Cauchy-Boltzmann theory or first-gradient theory. First, the 
basic kinematical relations and deformation measures are introduced. Then, the concept 
of internal contact actions, i.e. stress and heat flux are presented. Finally, the physical 
balance equations are outlined, along with their associated local forms, i.e. the Eulerian 
equations of motion. A more extensive and systematical presentation of these arguments 
can be found in e.g. |23, 61, 110, 17]. 


2.2.1 Kinematics of finite deformation 
Motion 


A continuum body is a set of material points that are in bijective correspondence, at each 
instant of time, with the geometrical points of a region of the Euclidean space, denoted as 
R?. The abstraction of continuum body is an approximation since the actual structure of 
the materials is discontinuous, but it is necessary for performing the classical mathematical 
operation (e.g. differentiation) and widely validated through experiments. 


Under the infinite possible configurations of the body, we call reference configuration Bo 
and current configuration B the configurations of the body at time t = 0 and t € Rt, 
respectively. The two configurations are equivalently called Lagrangian and Eulerian 
configuration and the material points are labeled as X and a, respectively. The generic 
nonlinear transformation that maps the reference configuration into the current one de- 
fined as 

x= (X,t) (2.5) 


is called placement. 


Deformation Gradient 


The second-order tensor 
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da = cof(F) dA 


Figure 2.1: Motion of a material body in the reference configuration Bo with position X 
and current configuration B with position x at time t. 


is called the deformation gradient. In the three-dimensional configuration, this tensor has 
nine independent components and characterizes the motion in the neighbor of a material 
point. 


The determinant of the deformation gradient characterizes the measure of the volume 
change. To exclude self-penetration, the Jacobian determinant must always assume a 
value greater than zero 


J = det(F) = AF xF):F>0. (2.7) 


With the definition of the deformation gradient and its determinant, it is possible to 
define the deformation of a surface and volume element from the reference to the current 
configuration with 


da = JFdA= cof(F)dA = SF x F)dA and dv=det(F)dVv. (2.8) 


The Jacobian determinant J is thus a measure for the relative change in volume during a 
deformation. Note that the tensor cross product operation, see |29], was utilized for the 
first time in the domain of solid mechanics in Bonet et al. |16]. 


The deformation gradient can be decomposed in a multiplicative form through the polar 
decomposition, and one can obtain: 


F=R-U=V-R (2.9) 


which is called right and left polar decomposition. R is an orthogonal tensor (i.e. 
RT = R!anddet(R) = 1) and it represents a rotation, and U and V are symmetric 
and definite positive tensors. These two tensors are called right and left stretch tensor, 
respectively, and are linked via the rotation tensor through the formula 


V=R-U.R. (2.10) 
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Deformation Measures 


Having defined the deformation gradient F and keeping in mind that it is not a suitable 
quantity for strain measures, we introduce the most common strain tensors related to 
either the material or the spatial description used in the field of nonlinear solid continuum 
mechanics. The two symmetric and positive definite strain tensors 


C =U? = FIF 


2.11 
b= R = FF" en) 


are known as the right Cauchy-Green tensor C, and the left Cauchy-Green tensor b, 
present in the reference configuration and in the current configuration, respectively. 


2.2.2 Stress tensors and heat flux 


In this section, we recall the concept of the internal actions, i.e. contact stresses and heat 
fluxes, based on Euler’s cut principle. Therefore we consider an arbitrary subdomain T 
cut out of the body B as illustrated in Figure 2.2. 


Figure 2.2: Euler’s cut principle, with the mechanical traction vector t = on, and the 
heat flux h=q-n 


Stress tensors 


The Cauchy traction vector t, acting within a given spatial point æ depends on the 
orientation of the unit vector n. The Cauchy theorem implies that a linear mapping 


exists 
t=on, (2.12) 


with the Cauchy stress tensor ø, which can be proved to be symmetrical (by imposing 
the balance of moments). The Cauchy theorem is also valid for the first Piola-Kirchhoff 


traction vector T, leading to 
T=PN, (2.13) 
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with the first Piola-Kirchhoff stress tensor P which in spite of the Cauchy stress tensor, 
is not symmetric. Another stress tensor can be defined through the Cauchy stress tensor 
o and the Jacobian determinant J 


T=Jo, (2.14) 
in which the tensor T is referred as Kirchhoff stress tensor. The second Piola-Kirchoff 


stress tensor S is defined as the pulled-back of the Kirchoff stress tensor 7. This tensor 
is a symmetric tensor that does not have a physical interpretation. It is defined as 


S = Pore: (2.15) 
The following relationships between the stress tensors defined above can be derived 


P = Jo FT 
S=F"P, (2.16) 
S=IPIGE=". 


Heat flux 


The Stokes heat flux theorem is the thermodynamic equivalent of Cauchy’s stress theorem 
(2.12). It postulates the linear mapping for the heat flux vector h on T 


h=q(a,t)-n, (2.17) 


with the Cauchy heat flux vector q, see Figure 2.2. The Lagrangian counterparts can be 
defined as 


H=Q-N with Q=JF'g. (2.18) 


2.2.3 Balance laws 


Balance laws are used to refer to the overarching principles or rules of nature. They 
apply to all types of materials and possess an axiomatic nature. They can be constructed 
in either a global integral form, affecting the whole body, or a differential local form, 
affecting each location on the body. Balance rules are used to relate an object’s external 
loads to its internal quantities. If these quantities remain constant during a process, they 
are referred to as conservation laws. The global and local balance principles of mass, 
linear momentum, angular momentum, and energy are briefly introduced in this section, 
as well as the entropy inequality using the first and second laws of thermodynamics. 
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Balance of mass 


With the spatial density field p(x, t) and the material density field pọ( X) the conservation 
of mass of an infinitesimal volume reads 


m= | wav = f pdv = const. (2.19) 


Bo B 


Using (2.8) the local form can be deduced as 


po = pd. (2.20) 


Balance of linear momentum 


The conservation of linear momentum is a generalization of Newton’s second law applied 
to a continuum body. As such, it is also known as the first Eulerian equation of motion 
and reads 


- [veda fbdo+ ftaa (2.21) 
B 


B oB 


Here, v is the spatial velocity field, b is an external body force field per unit volume (e.g. 
gravity). Using the Cauchy stress theorem and the divergence theorem, the local form is 
given by 


V-P+b-pw=0 (2.22) 


Balance of angular momentum 


The balance of angular momentum is the rotational analog of the linear momentum and 
is obtained by differentiating the angular momentum with respect to time 


d 
al exopdo= fe xbdo+ f æxtda, (2.23) 
B B aB 


leading to the local form 


o=0”. (2.24) 
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Balance of energy - The first law of thermodynamics 


The balance of total energy reads 


5 Goto udu = [(b-v+R)dv+ f(t-v- Ada, (2.25) 


B oB oB 


with the stored internal energy per unit volume u and a given heat source R. 


The local form 


d E 

ut FF V-Q+R (2.26) 
can be interpreted as the first law of thermodynamics, stating that the rate of change 
of total energy (kinetic and internal) of a system equals the total energy input into that 
system, conducted by mechanical and thermal energy. 


Dissipation postulate - The second law of thermodynamics 


The second law of thermodynamics is founded on the notion of entropy and describes 
energy changes in terms of their direction. For example, it defines that heat is always 
transported from a warmer to a cooler location within a continuum. In the subject 
of irreversible continuum mechanics, such as damage mechanics, entropy is a critical 
quantity. 


The second law of thermodynamics can be defined as 


T=8-Q>0, (2.27) 


with the total production of entropy T per unit time, the rate of change of entropy S and 
the rate of entropy input Q. It can be derived into the local dissipation postulate 


Din =P: F-96-v Q:v0>0, (2.28) 


known as the Clausius-Duhem inequality, with the free Helmholtz energy Y, the entropy 
n and the absolute temperature field 0. 


2.3 Shell kinematics 


This section presents the basic kinematics of Kirchhoff-Love shells used for the formulation 
of fiber-reinforced sheets in Chapter 5. Shells, defined as planar structural elements with 
a small thickness relative to their planar dimensions, are used in a variety of engineering 
applications. Plate and shell theories in continuum mechanics exploit this length scale 
disparity to reduce the three-dimensional problem of solid mechanics to a two-dimensional 
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one. The two most common formulations are the Mindlin-Reissner and the Kirchhoff-Love 
models. In the first case, the kinematics is determined by the displacements and rotations 
of the mean surface, while in the second case only the displacements are needed. With 
the advent of the isogeometric concept, see Section 2.5, the continuity required by the 
weak formulation of Kirchhoff-Love shells can be easily achieved. 


63 


Figure 2.3: Kinematics of the continuum shell in the undeformed and deformed config- 
uration. 


To be specific, we consider a three-dimensional continuum shell element as depicted in 
Figure 2.3, whose points in the Euclidean three-dimensional space R? are addressed by 
the convected coordinates (61,67, 6°). The reference placement X : By — R? is the 
embedding* 

X (0) =F(0%) +N (0%) , (2.29) 


where F := P(0%) is the parametrization of the reference mid-surface My C R°, which 
comes along with the covariant basis A, and the unit normal vector field N defined by 


Aa: =P a and N:= By Gs 


re a, 2.30 
| [A x Ad] ey) 


where Fa := 2£. Note that common operations defined on R? are extended for vector and 
90 P 


tensor fields in a pointwise fashion. For instance, the scalar field || Aı x Aa|| is determined 
by ||Aı x Ao||(@%) := ||(A1 (0%) x A2(0%)||, where we applied the cross product x of R? 
together with the standard Euclidean norm ||v|| := y/(v - v), v € R? induced by the inner 
product - of R*. In the reference placement (2.29), the coordinates are chosen such that 
0% € [—-h/2, h/2] denotes the coordinate in thickness direction with h as the thickness of 
the body. 


* Greek and Latin indices range in the sets {1,2} and {1,2,3}, respectively. Moreover, we will apply 
Einstein summation convention for repeated upper and lower indices. 
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The deformation of the body is described by the function æ : B — R?, which incorporates 
the Kirchhoff hypothesis by the kinematical ansatz 


(0) = r(6°) + 0’ n(0°) . (2.31) 


Here, r := r(0%) is the parametrization of the deformed mid-surface Q C R?, which comes 
along with the covariant basis ag and the unit normal vector field n defined by 


Aq =ra and n:= ———. (2.32) 


The kinematical restriction (2.31) allows to describe the body’s deformation exclusively 
by the geometry of the mid-surface Q. Therefore, we introduce some important geometric 
objects, that play a central role in the development of the mechanical theory. For a 
concise treatment of surface geometry, we refer to |24]. The co- and contravariant metric 
coefficients of the deformed mid-surface (2 are defined by 


das :=0G,.-ag and [a%?(6%)] := [aas], (2.33) 


where the square brackets indicate the component matrix. Consequently, it holds that 
aaga” = 64. The contravariant basis a” is introduced by the relation a” - ag = 03 
and can therefore be expressed as a? = a“? ag, which in turn leads to a°? = a®- af. By 
definition of the cross-product, we can rewrite ||aı x az ||? = [(a,-a1) (a2-a2)— (a1 -a2)°?] = 
det(aag). The partial derivatives of the covariant base vectors and the normal vectors can 
be expressed using the Gauss and Weingarten equations 


Tag = Qag = Uég Oo Fbagn and no = —dag a’, (2.34) 


where we use the Christoffel symbols of the second kind P'¿¿ and the covariant components 
of the curvature tensor bag defined by 


ag =A -Aap and bog :=@o,g°N. (2.35) 
Note that (2.34), can also be expressed as 
aag = Tas a? +dagn, (2.36) 
where 
Dago = Go daB = a?. af Aro = Toa Ado (2.37) 


denotes the Christoffel symbols of the first kind. Analogously, we can define the above 
introduced objects and relations also for the reference mid-surface Qo using the co- and 
contravariant metric coefficients Aag and A“P, the contravariant base vectors A“, the 
Christoffel symbols of the second kind reg as well as the covariant components of the 
curvature tensor Bag. 


The surface-deformation gradient F’, defined by dr = F dr, can be computed by inserting 
d9% = dO'A; - A% = 7,d0°- A% = dr - A“ into the expression dr = r,, d0% = 
aal dr - A?) = (a, € AC) dF and is consequently given by 


F=a,0 A". (2.38) 
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The covariant base vectors g; = Œ; of the three-dimensional body in the deformed con- 
figuration are 


Jal’) = (0%) +P n..(0°) and gu(0') = n(6") (2.39) 
and G; = X ; in the reference configuration 
Gal’) = A,(0%) + 6° Na(0%) and G3(6') = N(6°). (2.40) 


The covariant metric coefficients gi; = 9:9; of the three-dimensional body in the deformed 
configuration are determined as 


Gas(6') = 930.(0') =0 ’ (2.41) 


where we applied the definitions (2.33) and (2.34), together with the properties n-n = 1, 
a.:°n=O0andn,~-n = 0. Analogously, we obtain the covariant metric coefficients 
Gi; = G; : Gj in the reference configuration 


Gap (6") = Aag(0%) — 26° Bag(0%) + (6°)? Na(0%) - Nis(8*) , 
Gaz(0%) = Goal P) = 0, (2.42) 


G33(0) = 1. 
Neglecting the quadratic term in equations (2.41) and (2.42) yields 
Jap(9") = dag(0%) — 26% bag(O%) and Gaggl) = Aas(9*) — 29° Bag(0%) (2.43) 


assuming a linear behavior through the thickness of the body. Due to the block structure 
of the metric coefficient matrices, the contravariant base vectors gê and G" in the deformed 
and the reference configuration can be computed as 


g% =g% g, $=g=n, [I] = [ga], (2.44) 
and 
CHG" Gs, @ =G =N, EN =; (2.45) 


respectively. Similar to the surface-deformation gradient F', the deformation gradient of 
the three-dimensional body F is 
F= 4,8G', (2.46) 


which is easily obtained by de =xw.¿d0' = g(dX - GĦ) = (gi 9 GÌ) AX = F dX. 


the right Cauchy-Green deformation tensor for the surface is 


C:= FF = Ca A Q A’ , (2.47) 
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where we have defined Cyg(dag(0%)) := dag(0%). Introducing the difference 
Sie = Poe mi Do and SaBA I= Sas UA = Tapa = IE Usa 5 (2.48) 


it is convenient to define the second covariant derivative of the deformation with respect 
to the reference configuration as 


lab = Tag — ie a, = (Dis = IZ) as + bagn 


Bi o (2.49) 
= Sag Ac + bap N = Sapo a + bagn . 
Thus, we obtain 
(Tiag — Bag n) 9 A* Q Af =a, I I -NOK (2.50) 
with 
S7 := Si, A* O A? and K:= Kap A® Q A? = (Bas — bag) A* Q AP, (2.51) 


where in-plane and out-of-plane curvatures can be addressed by eight and four compo- 
nents, respectively. 


For the matrix material, the right Cauchy-Green tensor reads 


whose covariant components coincide with the metric components of (2.41); note that 
higher-order terms in 0% are neglected. Due to the dependence of the metric coefficients 
(2.43) on the first and second fundamental form, it is convenient to introduce the com- 
ponents of the right Cauchy-Green tensor as Ci; (dep (8%), Kap (8%), 0°) = gi¡(0%). The 
kinematical restriction in (2.31) does not allow for any change in thickness direction. For 
a fixed 0° = 63, the surface elements of the embedded surfaces in the reference and the 
deformed configurations X (6°, 6?) and x(0%, 6°), respectively, are 


da(0%, 0%) = ||g1 (0°, È) x 92(0°, 03) | dO! dé? = \/det(gus(0°, 03) do! de? , (2.53) 
dA(6%, 6°) = ||G (0%, 6°) x Ga(0%, 6°) || do! dé? = 4/det(Gag(0%, 03) dé d0? . (2.54) 
Consequently, the in-plane Jacobian determinant 


det ( det (Jag (0%, 07) (8°, 63)) 


det (Gas(0°, 03) a 


can be interpreted as the areal dilatation, i.e. the ratio between the deformed and the 
reference area elements da and dA, respectively. 


2.4 Phase field approach to fracture 


The modeling of fracture phenomena with phase-field methods have gained increasing 
attention in recent years. In this section we will introduce the underlying concepts for the 
case of non-linear brittle fracture. In subsequent chapters an extended hybrid model will 
be deduced for porous-ductile fracture in fiber-reinforced thermo-elastic-plastic solids. 
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2.4.1 Crack regularization 


Bo OBo 


Figure 2.4: A schematic representation of an approximation of an internal discontinuity 
by a phase-field 


We consider an arbitrary elastic body of the bounded domain By in the reference con- 
figuration and an evolving internal discontinuity boundary I , that represents a set of 
discrete cracks, as shown in Figure 2.4. In accordance with the classical brittle fracture 
approach of Griffith [55], crack nucleation and propagation arises upon the attainment of 
a critical local fracture energy density ge. The total potential energy of the body then 
reads 


wet = W* + Wf = [wav fs. dA, (2.56) 
Bo T 
whose minimization can describe the entire fracture process, controlled by a potential 


elastic strain energy term and a term related to the work required for the creation of new 
surfaces. 


The tracking of the crack discontinuity boundaries in classical fracture approaches often 
requires complex and costly computations, particularly when considering complex shaped 
cracks or even crack branching. The underlying concept of the phase-field formulation is 
to approximate the discrete fracture surface with a continuous scalar field, the so called 
phase-field s( X,t), where the value s = 0 refers to the undamaged and s = 1 to the fully 
broken state of the material. Thus, the crack surface energy can be approximated as 


fdas fudsyav = [arsav. (2.57) 


r Bo Bo 
with a crack density functional y, typically defined as 


ee 5 V(s) -Vís). (2.58) 


The smooth bases of isogeometric analysis presented in Section 2.5.2 provide a convenient 
computational framework to apply higher-order regularization and thus obtain better ac- 
curacy and convergence rates, see Borden et al. [19]. The common fourth-order approach 


reads 
qls, V(s), A(s)) = a” + Lvs) -V(s) + E Als) A(s). (2.59) 
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The length scale parameter / in the the crack density functionals characterizes the width 
of the smeared crack profile. For s 0, the approximation converges to the sharp crack 
T. It can also be considered as a material property governing the critical value of stress 
required for the crack initiation or growth. For a one-dimensional bar under tension, | 
can be derived with the critical tensile stress o., as 


_ 27E Gg. 

256 0.2 
see |13] for details. Though an extension in higher dimensions and complex mixed-mode 
fracture analysis is difficult, this can give a rough estimation on how to choose l. In 


order to correctly resolve the steep gradient of the smooth crack profile, a fine element 
discretization should be employed in the region where a crack is expected. 


(2.60) 


---- 2Morder 
4% order 


Ez 
l 


Figure 2.5: Analytical solution for unidimensional phase-field profiles for second- and 
fourth-order formulations 


The unidimensional analytical solution of the corresponding Euler equations of the crack 
density functionals (2.58) and (2.59), after minimizing So, ge y(s) dV and imposing the 
initial condition s(0) = 0, are depicted in Figure 2.5. Note that the higher-order formu- 
lation has a slightly narrower profile, meaning that a smaller region needs to have a fine 
mesh in order to correctly resolve the crack. For this and for the abovementioned reasons, 
the fourth-order model is adopted for the simulations in this work. 


2.4.2 Elastic strain energy 


Next, the loss of material stiffness in the crack is incorporated via the degradation of the 
strain energy through the phase-field. The simplest constitutive description would be 


W°(F,s) = g(s) V(F) (2.61) 


with the linear degradation function g(s) = (1 — s$). Assuming that fracture requires 
a local state of tension, a physically more reasonable formulation distinguishes between 
crack-insensitive compressive (~) and crack-sensitive tensile components (+). A common 
additive decomposition of the strain energy density reads 


W°(F,s) = g(s) V*(F) + Y (F), (2.62) 
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where the degradation induced by the phase-field acts only on the tension term of the 
strain energy function. A possible formulation for large deformation, proposed by Amor 
et al. [8] and Borden |18], is based on a the following decomposition of the strain energy 
into volumetric and deviatoric parts 


ped = a 
VU = Pval J) 


w= WVaev(F) F Pval J) 
yv” =0 


(2.63) 
elseif J>1 


Multiplicative decomposition of crack-sensitive and -insensitive components have also 
been proposed over the principal stretches A, for the non-linear regime. Hesch et al. |126] 
introduced an anisotropic split via A, = Az At with 

(A. — DI 


ELA. —1 
At = O Deon +1 (2.64) 


leading to the fracture insensitive part of the deformation gradient 


F= Y (at) A na @ Nay (2.65) 


a 


where n, and N, denote the principal directions of the left and right stretch tensors, 
respectively. Since the just mentioned decomposition is not compatible with a standard 
elastoplasticity formulation, which relies on different mechanisms for the deviatoric and 
volumetric contributions, we introduce a novel multiplicative decomposition in this work, 
defining the fracture insensitive, isochoric part of the deformation gradient and the frac- 
ture insensitive of the Jacobian determinant as 


= - JJ® if J>1 
F=) (JA On 0N, and J= ; 2.66 
2 ue an J else ( ) 
The associated elastic strain energy reads 
W°(F, 5) = Veo (F(F,5) + Val J(F,8)) (2.67) 


for which typical decoupled strain energies can be taken into account. 


Note that more alternative formulations for the tension-compression split have also been 
presented, like decompositions based on the symmetric and antisymmetric parts of the 
strain tensor |51] or with respect to the crack orientation |119]. 


The degradation function g(s) controls the transition of the behavior of the material from 
the intact to the cracked state and has generally the following properties 


g(0)=1, g(1)=0 and g(1)=0. (2.68) 
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Commonly, the following quadratic degradation function is considered 
g(s) = (1-5) +e. (2.69) 


The positive small factor e = 0 can be used for avoiding zero stiffness of the material in 
a fully cracked state. With a more general cubic degradation function in the form 


g(s) = (ag — 2) (1 — 5)? + (3 — ag) (1-5)? + € (2.70) 


a small but positive slope of the function at s = 0 allowing better crack nucleation 
is achieved with small positive values for the degradation modeling parameter ag. For 
brittle fracture, a perfectly linear stress-strain relationship, i.e. a pure elastic behavior 
with no softening prior to crack nucleation can be observed with with ag —> 0, see |18]. 


2.4.3 Balance equation of the phase-field 
By varying the density functions of the stored energy with respect to s, we obtain the 


additional constitutive quantities 


H = —6,0° = -9, 0° 
rh =5,0! =0,0! — Div[0y.V/] = Es — ge lp As, 


namely the driving force of the phase-field H and the crack resistance force rf. Note 
that the second order functional I is considered here. The strong form of the mechanical 
system can then be extended with 


H= £s — gely As. (2.71) 


2.5 Isogeometric analysis 


In this section we provide an introduction of the isogeometric concept used for all geomet- 
rical discretizations in this work. In conventional finite element analysis, the spline-based 
CAD geometry is approximated with low-order, usually linear, Lagrange polynomials ba- 
sis functions. The term "isogeometric analysis" implies that the FEA model employs the 
same spline-based description as in the CAD geometry. Additionally, it provides higher 
continuity of the finite element approximation, needed for certain models, i.e. the fourth- 
order phase-field model, the Kirchhoff-Love theory and the strain-gradient formulation of 
the fiber material. 


20 2 Fundamentals 


2.5.1 B-splines and NURBS 
Univariate B-splines 


Univariate B-spline functions are piecewise polynomial functions with compact support. 
They are defined in parametric space using a so-called knot vector denoted = . In a one- 
dimensional space the knot vector is a set of finite and monotonically increasing sequence 
of real numbers = = [£1,€2,...,En+p+1), where p is the polynomial degree and n is the 
number of basis functions of the B-spline. 


Using the convention that a fraction in front of the basis functions is set equal to zero in 


the case of the denominator being zero (i.e. 2 := 0), the B-spline functions are defined 
recursively on the vector of knots with the Cox-de Boor formula [27] as 
Base ¿BOE Pl: (2.72) 


= = 
Cap E + &rprı — S41 


beginning with 
; 1 if <E< Gist 
Bolg) = (2.73) 


0 otherwise 


for p = 0. It can be noted that for p = 0 and p = 1, the basis functions of isogeometric 
analysis are identical to those of the standard piecewise constant and linear finite elements, 
respectively. 


Each knot represents a coordinate value in the parametric space. If the knot vector is 
chosen equal to a set of following integers, we refer to natural B-spline. If the knots are 
distributed uniformly, the vector of knot is said to be uniform. It is said to be open 
if the first and last knots are repeated p + 1times. Open knot vectors are standard in 
the CAD literature. In one dimension, basis functions formed from open knot vectors 
are interpolatory at the ends of the parametric space interval, [€1, Emyp+1), and at the 
corners of patches in multiple dimension but they are not, in general, interpolatory at 
interior knots. An example of quadratic basis functions for an open, uniform knot vector 
is presented in Figure 2.6. 


In general, basis functions of degree p have p — 1 continuous derivatives. If a knot is 
repeated k times (with k < p), then the number of continuous derivatives decreases by 
k. When the multiplicity of a knot is exactly p (i.e. k = p), the B-spline basis functions 
are C° continuous at that knot point. This is illustrated in Figure 2.6 at the location 
of the repeated knot € = 4. Furthermore, if the multiplicity of the knot point is p+ 1 
(i.e. k = p + 1) the B-spline basis functions are discontinuous at that knot point and as 
a result two separate B-splines patches are formed. 


In the perspective of finite element analysis, the univariate B-spline basis functions, as 
well as the multivariate B-spline and NURBS basis functions have important features, 
namely 
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< 09 


0 
0 1 2 3 4 5 


Figure 2.6: Quadratic basis functions for the open-uniform knot vector 
€ = [0, 0,0, 1, 2,3, 4, 4, 5,5,5] with varying degrees of regularity. 


e Partition of unity > Bi(é)=1 


e Nonnegativity over the entire domain BE )>0 
e Compact support of each Bi (£) in [E €i+p+1), 


e Linear independence > B¡(€)c,=0 & c=0 


The derivatives of B-spline functions are represented in terms of lower-order B-spline 
basis. Thus, the first derivative of the ith B-spline basis function of degree p is given by 


dp oP 
© pile) = 
ag Ph) = Be 


and consequently the jth derivative with 


di 7 = oP —_ ai! i p di a4 
jer PHO = Pe (ger 0) ae reno) em 


Multivariate B-splines 


By 1(8) — Be (2.74) 


One can easily extend the previous concepts to define multivariate B-splines for dimension 
d € {2,3} by using the dyadic product E = =; 08... 8 E4 of univariate knot vectors 5,, 
with the direction in the parameter space i € {1,...,d}. More precisely multivariate 


B-splines are defined as 
d 


Bo = Bie) =] | Bele), (2.76) 


l=1 
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with the associated univariate B-spline basis functions Bn (E!) where the multi indexes 
i = [iz,...,%q] and p = [p),..., pa] represent the position in the tensor product and the 
polynomial degree in the parameter space, respectively. The assignment of the global 
numbering index G to the multi index ¿ can be found [26]. 


Using (2.74) and (2.75), the jth partial derivative of the Gth multivariate B-spline with 
respect to £! can be obtained as 


of BS ok 
Da zen B ojl Bi (é!), (2.77) 


A 


NURBS 


B-splines are convenient for free-form modeling, but they lack the ability to exactly repre- 
sent some simple shapes such as circles, cylinders and ellipsoids. This limitations can be 
overcome with the extension to non-uniform rational B-splines (NURBS), which enable 
more control over local geometry shape due to non-uniform treatment of each control 
point influence on the curve. This is obtained by the use of weights w; for each control 
point respectively according to 


d 
TB ls )w 
G_ yie _ = 
NS = Ni(é) = = eaves (2.78) 
2, H BACE) w; 
Jai 
Note that by setting the same value for each weight we obtain the original B-spline. 
By applying the quotient rule, the partial derivative of N“ with respect to €! reads 
BB we OBA _ BA OBB wp 
one % ar Bar 


an” 2 
(Gem) 
B 


Higher jth order derivatives can be obtained by 


ði G gr=m G om B 
gel = Y BB wg i (2.80) 
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Hierarchical refinement 


The univariate and multivariate B-splines and NURBS can be refined by different ways 
without changing the shape of the geometry. 


Through knot insertion, the solution space is enriched by adding more basis functions 
of the same order while keeping the spline unchanged. Insertion of new knot values clearly 
has similarities with the classical h-refinement strategy in FEA. 


In the process of order elevation, the B-Spline/NURBS-order is raised by increasing the 
multiplicity of each knot (and without adding any new knots values). It can be thought 
of as p-refinement in FEA. 


Lastly, hierarchical refinement of multivariate B-spline and NURBS bases have been 
proposed in |35] with preservation of the global properties in terms of smoothness and 
continuity of the unrefined mesh. 


2.5.2 Geometric discretization 


In the isogeometric analysis framework, the computational domain is not approximated, 
but exactly described. To be specific, we introduce a non-linear transformation from the 
parametric space B into the physical space B 


T: BoB 


Gr (2.81) 


By associating the NURBS basis functions NC with a net of control points Xg € R“ 
related to the geometry under consideration, i.e. a curve or surface as seen in Figure 2.7, 
we obtain 


G=T(E) = Y N Xe =X Ni(E)X,. (2.82) 
G i 


Figure 2.7: Physical mesh and control mesh of a quadratic NURBS curve (left) and a 
quadratic NURBS surface (right). The red points denote the control points. 
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Accordingly, rational approximations of the deformed geometry y and its variation dp 
are defined as 
p =X N Xe, bp = Y N%Xo, (2.83) 
G G 
which we can insert into the corresponding terms of the balance equations to obtain the 
discrete formulation of our problem. 


3 Phase-Field Modeling of Non-Linear 
Thermo-Porous Ductile Fracture 


In this chapter, we present a framework for simulating porous-ductile fracture in isotropic 
thermo-elastic-plastic solids subjected to large deformations. The proposed model is 
based on a triple multiplicative decomposition of the deformation gradient and an expo- 
nential update scheme for the plastic return map in the time discrete setting. In addi- 
tion, a modified Gurson—Tvergaard—Needleman GTN-type strain-gradient thermoplastic- 
ity model is combined with a second-gradient phase-field fracture approach to account for 
a temperature-dependent growth of voids on microscale followed by crack initiation and 
propagation on macroscale. The multi-physical formulation is completed by incorporating 
an energy transfer into the thermal field such that the temperature distribution depends 
on the evolution of the plastic strain and the crack phase-field. Eventually, this physi- 
cally comprehensive fracture formulation is validated by experimental data and applied to 
complex three-dimensional geometries. Based on the available comparative experimen- 
tal data, the specific material models and parameters in this chapter are provisionally 
adapted for steel-like materials. Generally applicable to thermo-elastic-plastic solids, the 
framework and formulations presented here are specifically adapted to the final polymeric 
matrix material in the subsequent chapters. 


3.1 Configuration and kinematics 


As introduced in Chapter 2, we consider a body of interest By € R? in its reference 
configuration, along with the deformation mapping (X,t), the material deformation 


gradient F(X, t), the crack phase-field s(X,t) and the absolute temperature 0(X, t). 


Regarding the plasticity, we introduce the equivalent plastic strain a(X,t) and its dual 
dissipative resistance force r?(X,t) as 


a(X,t):Byx TOR and r?(X,t): Box TOR. (3.1) 
The plastic deformation gradient F? represents an additional strain measure field 


F°(X,t): Box TR, with JP =det[F?] > 1. (3.2) 


Within the framework of a GTN-type plasticity model, the void volume fraction f is 
introduced as a micromechanically motivated damage parameter for growing cavitation 
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in plastic solids 


volume of voids 


(X,t): Box T > [fo fe] with f= (3.3) 


total volume 


where fo is the initial porosity and fr is the final void volume fraction. The above 
introduced variables characterize a multi-field setting of porous-ductile fracture in thermo- 
elastic-plastic solids based on five independent fields 


U :=[p,5,0,0,rP), (3.4) 


the finite deformation map q, the absolute temperature 0, the crack phase-field s, the 
equivalent plastic strain a and the dissipative plastic resistance force rP. Note that the 
Lagrangian plastic deformation map F? and the void volume fraction f will be condensed 
within the balance equations. 


As is common in nonlinear elastoplasticity, we apply a multiplicative split of the defor- 
mation gradient and its determinant, and obtain the elastic components 


F?=F(FP? and P=J(P). (3:5) 


Postulating that fracture requires a local state of tensile/shear deformation as discussed 
in |8, 88, 37, 40], we define the fracture insensitive, isochoric part of the deformation 
gradient and the fracture insensitive of the Jacobian determinant as follows 


MERO if Je>1 


Je else 


Fe = ye n, ON, and = (3.6) 
The degradation function g(s) is often defined by the following cubic function g(s) = 
(ag — 2) (1—5)? + (3 — ag) (1—5)?, where a, € (0, 2] is a degradation modeling parameter. 
Note that we use a simplified quadratic degradation function (1 — s)?, which is obtained 


by setting ag = 2. Moreover, AS = (A8)9®) are distortional stretches, defined via the 
isochoric, elastic parts of the principal stretches A? = (J°)~/2)¢, and finally, the vectors 
Na and N, represent the principal directions of the left and right stretch tensors, see 


[37]. 


Furthermore, in order to subsequently derive objective (frame-independent) constitutive 
relations, we introduce the elastic left Cauchy-Green 


b° = F (CP)! FT (3.7) 


with CP = (FP)! FP. Consequently, the fracture insensitive part of the elastic left 
Cauchy-Green tensor reads 


bo = F°(F)? = YOA? na @ ma. (3.8) 


a 


In contrast to the isochoric von Mises plasticity using J = y/det[b*], the volumetric plastic 
void growth considered in this work leads to the total determinant 


J = \/det |b] y det[Cr). (3.9) 
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3.2 Variational formulation 


Next, we propose the constitutive framework for the thermomechanical damage model. 
To be specific, we introduce constitutive energetic and dissipative functions and derive the 
necessary relations and evolution laws to formulate the multi-field variational problem. 


3.2.1 Energetic response function 


By coupling gradient thermoplasticity with gradient damage mechanics, the density func- 
tion of the stored energy W for the multi-field problem is defined by 


T = (F, F”, 5,9) + 0°(#) + UW (a, Va, 0) + Vs, Vs, As). (3.10) 
The elastic contribution is composed of volumetric and deviatoric parts as 
Te = VETAS, JP, 5), 0) + VO, (FE, PP, J,s)) (3.11) 


where, without prejudice to the generality, we consider in the examples an extended 
Neo-Hookean material law given as follows 


u =5 (== - 1029) -30n0 0) (FF) 638 


and se z en, 
De. (F°) = E (F. : F° — 3). (3.13) 


Therein, 4 > 0 and « > 0 denote the shear modulus and the bulk modulus, respectively. 
00 is a reference temperature and 3 is the linear thermal expansion coefficient. Moreover, 
for a compact representation we introduce the abbreviation F° : F° = Y (A? )?. Next, the 


a 
purely thermal contribution to the energy (3.10) is assumed in the simple form 


v°(0) = (00 — din ($) (3.14) 


where c > 0 is a constant parameter representing the specific heat capacity. Satisfying 
the heat conduction inequality via Duhamel ’s law of heat conduction, we define the Piola- 
Kirchhoff heat flux vector as 


Q(F,s,0,V0) := —K(F,s,0)V0. (3.15) 


Phenomenologically motivated, we formulate the material thermal conductivity tensor K 
in terms of the phase-field parameter s 


K(F,5,9) := [Ko (1 — wx (8 — %)) 1—s) + KO s] C7}, (3.16) 
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so that when a fracture occurs, the conduction degenerates locally and leads to a pure 
convection problem where the heat transfer depends on the crack opening width. More- 
over, wk is a thermal softening parameter, Ko is a conductivity parameter related to the 
reference temperature, K°®Y is a convection parameter and C = FT F denotes the right 
Cauchy-Green tensor. 


Based on the particular formulation for gradient plasticity in |2, 5, 37, 35], the plastic 
energy contribution is defined in terms of the equivalent plastic strain a, its gradient Va 
and the temperature field @ as 


> f 2 
D(a, Y 0,6) = foan då + (0) È Val. (3.17) 
0 


As described in |2], lp is a plastic length scale associated with a strain-gradient hardening 


effect that takes size effects into account in order to resolve the nonphysical mesh sensi- 
tivity of localized plastic deformation in softening materials. The temperature-dependent 
isotropic local hardening function y(a, 9) can be determined experimentally for the re- 
spective material. In the context of metal plasticity we use this particular saturation-type 
function outlined in [112] 


(a, 0) = Yoo(#) — (Yoo(#) — yo(0)) exp[—wpa] + h(9) a, (3.18) 
with the three temperature-dependent material parameters yo > 0, Yoo > Yo and h > 0 


defined as 
Yyo(9) = yo() (1 — wo(9 — 00), 
h(0) = (00) (1 — wa (0 — 00)), (3.19) 
Yoo(9) = Yoo) (1 — (0 — 00). 
Here, wp is the hardening/softening parameter, wo is the flow stress softening parameter 
and wp is the saturation parameter. The initial yield stress yo determines the threshold of 
the effective elastic response. The variational derivative of WP with respect to a yields 


rP := ôa DP = ô, — Div[Oy, DP] (3.20) 
reflecting the characteristics of the gradient-extended model under consideration. 
Lastly, following the phase-field approach to fracture introduced in Section 2.4.3 , the 
phase-field fracture contribution is given in terms of the crack-density function as 


T! (s, Vs, As) = Ge(a) A(s, Vs, As) 


Fela) 2, Fla) le Fla) lE aa (3.21) 
= == —— Vs. Vs + ——— (AsV. 
11, 5 + 5 s-Vs+ 1 (As) 
In the context of ductile fracture, we first require that lp > lẹ such that the regularized 
crack zone lies inside of the plastic zone. Secondly, we propose a degradation of the 
Griffith-type critical energy release rate Je. In Miehe et al. |83], YP (a, Va) is degraded as 
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well and thus, included in the driving force of the phase-field. This behavior can be equiv- 
alently described by decreasing the local fracture energy by a factor of 1 — me g (s) Ur. 
Physically interpreted we can state that the critical fracture energy is decreased by the 
amount of energy included by the plastic work density function. Regarding the exponen- 
tial character of WP, we therefore propose the following degradation of the critical fracture 
energy density 


Fela) = Ic,p ag Ic,e exp[—wral, (3.22) 
using the modeling parameters {9.,¢, Jc,p; Wt }- 
Finally, the crack resistance force can be obtained by the variational derivative with 


respect to s 
rf := 6,0! = 9,VÍ — Div[Oy,W*] + Alan, Y]. (3.23) 


3.2.2 Dissipative response function 


Regarding the porous elastoplastic material behavior, we consider a GTN type function 
based on |57, 123, 96], which implicitly defines the microscopic effective stress 5 := (T) 
in terms of the Kirchhoff stress r =0/J and the void volume fraction f 


C2 
Y“(r,5) := > + 2q1 fcosh Gal - (1+ (11?) =0, (3.24) 


[02 


Here, 0 = 3/2 ||Taev/J|| represents the von Mises equivalent stress, whereas p = 
str[r/ J] denotes the local pressure along with the growth-based void volume fraction 
f and fitting parameters qq = 1.5 and q2 = 1.0, see |6, 85]. Note that for qq = 0 
the influence of the pressure and the void volume fraction vanishes, i.e. in this case the 
Gurson-type yield criterion is identical with the classical von Mises yield criterion with 
7 = eq. Focusing on void growth and thereby neglecting other influences such as void 
nucleation or void softening due to shear, the evolution form of the void growth reads 
f = (1- f)tr[d?]. In line with |85], the current void volume fraction is given in terms of 
the plastic deformation as 

1-0 

JP 

with the initial void volume fraction fo. With the effective stress (7) and the dissipative 
resistance force rP we define the plastic yield function as 


f=1- (3.25) 


a 


lölr),r?)=öl(r)-r?. (3.26) 


A plastic Lagrange multiplier AP is introduced to enforce the Karush-Kuhn-Tucker con- 
ditions 
OP <0, APOP=0. (3.27) 
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For the incorporation of the fracture mechanical behavior, we define crack threshold 
function as 


D(H -r =H- rf, (3.28) 


where the energetic driving force H is bounded by the crack resistance force rf dual to the 
fracture phase-field s. Similar to plasticity, we introduce a fracture Lagrange multiplier 
AT to enforce the Karush-Kuhn-Tucker conditions 


M>0, 6 <0, MB =0. (3.29) 


3.2.3 Evolution equations 


The evolution of the thermo-elastic parts of the energetic function is given by 


Te To 
d= dra, a Ow awe (8 +09) 
Lgo. [ve ul = be i 6 
di aie Be pee rag ey) 
with the material time derivative 
b° = 1b° + bo + F(OP) FT, CP = (PAE? (F?)F. (3.31) 


In case of isotropy, the skew-symmetric part of the spatial velocity gradient vanishes, i.e. 
l = d, and O° commutes with b° such that the first term in (3.30) can be written as 


oye. jaw awe aw l 
b= [El a l] a b° | : [F(ĊP) -FT (b°) 
ðbe Obe 35 Obe + öbe [ (C ) ( ) ] 
ze (3.32) 
owe 
= |2 b| : [d — d?], 
| la 0, 
where d?” = -1F(OP)-!FT(b°)-! denotes the Eulerian plastic rate of deformation ten- 


sor. 


Regarding the partial derivatives in (3.30), we introduce first relations related to the 


Kirchhoff stress 
T = Tdev + Tool 


= OW bev he Owe i b°, (3.33) 
obe Ob* 
along with the driving force of the crack phase-field 
awe 
H = -—, 3.34 
Os ° oe) 
and the specific entropy 
a (t + vo) 


(3.35) 
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Based on the concept of maximum dissipation, we define an extended dissipation potential 
and obtain a constrained optimization problem as 


V:= sup sup |r: d? — J(1— f)r?a + (H -— r) — PDT, rP) — MOH — r’), 
x 


T,rP,H—rÍ AP,A 


(3.36) 
where the Lagrange parameters AP and Af control the non-smooth evolution of the plastic- 
ity and the fracture, respectively. This allows us to define the associated plastic evolution 
equations as follows 


age a» Oo 
d? = AP — d a=—-—~—, 3.37 
ar NS ET) Or’ en 
as well as the evolution equation for the crack phase-field as 
aot 
x 3.38 
s= OGL — rh (338) 


along with the loading-unloading conditions introduced in (3.27) and (3.29). A penalty 
regularization of the Lagrange multipliers can be utilized as follows 


P= cu) 30 ad wee (OH = e >0, (3.39) 


where 7, and nr are additional material parameters which characterize the viscosity of 
the plastic deformation and the crack propagation. Moreover, the Macaulay brackets 
are defined by (x) := (x + |x|)/2. Eventually, we define the internal dissipation density 
function Din: as follows 


Du := Up T : dP? + HS, (3.40) 
where v, is a constant dissipation factor typically chosen in the range of 85% to 95% in 
the context of thermoplasticity, see |112, 130, 75]. In addition, 1; is introduced as fracture 


dissipation factor based on the discussions related to an energy transfer into the thermal 
field in |38, 107] and the references therein. 


3.2.4 Strong and weak form 


The strong form of the coupled problem is summarized in Table 3.1. 


32 3 Phase-Field Modeling of Non-Linear Thermo-Porous Ductile Fracture 


1. Stress equilibrium 


div[o] +b = 0 3.41) 
2. Energy balance 
61) + Div[Q] — Din -R = 0 3.42) 
3. Fracture force as 
SU + (Hr!) =0 3.43) 
4. Hardening force pi 
ôaYP — rP = 0 3.44) 
5. Plastic force i 
owe 
-2 be =0 A 
ab +7 3.45) 
6. Plastic strain Es 
pp l a 
dP eee =0, AP = —(0P) 3.46) 
OT Np 
7. Equivalent strain 
AP O®P 
o eres Sy 3.47 
JOA oP ) 
8. Phase-field pA 
ao! les 
5-A =0, A of 3.48 
OH FF) me 


9. Dirichlet and Neumann conditions 


p = P(X, t) on OBS" and TFYN =T(X,t) on OB?" 


0=0(X,t) on aBee and - Q- N = Qy on IB 

2 (3.49) 
s=10n0B5* and Y (Zas -s) -N = 0 on 087" 
a = 0 on OBS" and Va- N = 0 on 08,” 


10. Initial conditions 


P(X,0) =o, $(X,0)= v0, 0(X,0)=00, s(X,0)=0, a(X,0)=0, r?(X,0) =0 
(3.50) 


Table 3.1: Strong formulation of the coupled problem. 


Moreover, the admissible test functions related to an extended set of global primary fields 
AL is given as 041 := [0(p, âs, 00, Oa, Or? } i.e. the variations of the deformation, the crack 
phase-field, the absolute temperature and the isotropic hardening along with its dual 
driving force. Their spaces are defined as 

V? = {dp € H'(Bo) | dy = 0on OBS}, 

V° = {5s € H?(B,)|ös =0 on dBi A Vös- N ondBo}, 

V? = {50 e H!(B,)|60 = 0 on ðB}, (3.51) 

V° = [da € H'(Bo)|da = 0 on Bgt}, 

) 


| 
vr = {5r? e £L?(Bo)}, 
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where H* and k > 0 denotes the Sobolev functional space of square integrable functions 
and derivatives. Note that the phase-field is required to be in H?(B,), which has substan- 
tial effects for the numerical solution of the coupled problem using finite elements. The 
weak form of the balance equations of the coupled phase-field approach to porous-ductile 
fracture in non-linear thermo-elasto-plastic solids reads 


Go = ir: Vibe) - 9: B| dV — if 6p: TdA=0, 
Bo aBy” 
a FB 
G; = [ös ne $ — ôS Xt (x- des) + Xt Ge le Vs. Vós + TEL AsAós| dV =0, 
Bo i 
Go = fw (Où — Din — R) — Q - V56] AV + i 50 H dA=0, 
Bo OB,” 
Ces i. [Sa (9-7?) +y(0) È Va- Véa] dV =0, 
Bo 


Gro i= for (mi ne) dV =0. 


Bo 


(3.52) 


Here, b and B are given body forces per unit volume of the deformed and reference con- 
figuration, respectively, N is the outward unit normal vector and T the surface traction 
at the boundary. Moreover, R represents the heat supply per unit volume and H = Q-N 
the heat supply on the thermal Neumann boundary. Appropriate Dirichlet conditions are 
given in Table 3.1, where (fp denotes a prescribed deformation and @ a prescribed temper- 
ature. Additionally, the Dirichlet boundary condition for the crack phase-field is given 
by 

s=1 on OB EN, (3.53) 


ensuring that a broken state remains broken. The Karush-Kuhn-Tucker conditions in 
(3.27) and (3.29) are evaluated by inserting local variables defined as 


1 for £30 1 for ®>0 
Xe =: and Xp =: ; (3.54) 
0 otherwise 0 otherwise 


Note that by using the local variable xr in (3.52)2, we demand $ > 0 for thermodynamical 
consistency, avoiding a transfer of energy from the fracture phase-field into the mechanical 
field. This prevents healing effects, which may be taken into account as well. We can 
also set xf = 1 and restrict only the fully broken state, i.e. we allow for healing until the 
phase-field reaches one. 
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3.3 Spatial discretization 


Following the isogeometric analysis approach presented in Section 2.5.2 and the approxi- 
mations of the deformed geometry ¢ and its variation dy defined there, the approxima- 
tions of the phase-field s and the temperature field 9 as well as their variations 6s and 60 
read 


=) Rósa, 05 =D) R654 (3.55) 
AET AET 
and 
=) Ra, 00" = Y R4804. (3.56) 
AET AET 


Therein, we use the same spline based approximations for all fields, satisfying the required 
continuity s",ds" € H?(Bo) of the fourth-order phase-field approach, where the global 
shape functions R^ : Bọ > R are associated with control points A € Z = {1,..., N} 
with the total number of control points Y. 


For the hardening variable a and its variation da, as well as for its dual driving force r? 
and its variation dr?, we make use of linear shape functions N“ defined on the physical 
mesh representation of the NURBS geometry with nodes a € J = {1, ..., n} with the 
corresponding number of control points n: 


a” =X N%o,, da? = X Nba, (3.57) 
aET aET 
and 
reh = XO Ner, rah = NO (3.58) 
aET aET 


The more natural choice of approximating the two plastic fields with the same quadratic 
NURBS shape functions leads to oscillations within both fields, indicating stability issues 
as already observed in [7]. The above described scheme has shown to be stable and 
numerically robust within our numerical examples. 


3.3 Spatial discretization 35 


Inserting the approximations into (3.52) yields the semi-discrete set of coupled equations 


Gt = 0q4: pore dV - Feta 


G? := 684 oa [rw dV+K2Bsg|, 


Bo 


G} := 604 | | (m RAR? 0, — RADE, — VRAQ*) dv - QA], (3.59) 


= ÓN 5 — Nr?) AV + Ko], 


Aph/-h „ph 
GY, := ôr? | May — J xpi? EA ay 


Gh 


a 


The Kirchhoff stress tensor, the local entropy and the phase-field driving force are given 
by 


Aven (Bem sh, 6h) 
Dy. 2 ame) por ¡PA 
J Ober on 


AW (Be, sh, 0%) 
Osh 


AVEO (Beh, sh 9») 
h ’ 
(3.60) 
respectively. The semi-discrete definition of the heat flux vector along with the conduc- 
tivity tensor read 


and yb = — 


Q? = K*VR404 and K® = [Ko (1—wx (0% —0p)) (1—5?) + Km gh] (C>), (3.61) 


where 


Ch = q4 - qg VR48 VRE. (3.62) 


Moreover, the semi-discrete external contributions in (3.59); and (3.59)3 are formulated 
as 


Pea — i, RAB dV + J RT dA and Q™*4 = / RAR AV + / RAQn dA. 


Bo OB" Bo OBJ" 
(3.63) 
The coefficients of the matrices in (3.59), are given by 
MP = m f RAR? dV, 
(3.64) 


1 1 1 
KM? =~ i IX (quer +BVRA.VRB4+ ARMAR”) av, 
f 


Bo 
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whereas the matrices in (3.59)4 and (3.59); take the form 


KË = yl j VN@.VN?dV and M? =n, / N®N? dV. (3.65) 


Bo Bo 


The semi-discrete definitions of the critical fracture energy density and the local hardening 
function read 


g = Ic,p F Ic,e exp[-wra®] (3.66) 


and 


J = Yoo(O") — (Yoo(O") — yo(6")) exp[—wpa*] +h(0%) a”, (3.67) 


respectively. The discrete void volume fraction optained from f? = 1—(1— fo)/JP* with 
JP» = ,/det[CP*] and J" the discrete version of the Jacobian determinant. Eventually, 
the internal dissipation in semi-discrete form reads 


Dey = Vp T” : dP” + up Hs, (3.68) 


int 


using the semi-discrete evaluation of the plastic deformation rate tensor d?”. 


3.4 Temporal discretization 


Finally, in order to obtain a set of non-linear algebraic equations to be solved via a Newton- 
Raphson method, we discretize the semi-discrete coupled problem (3.59) in time. There- 
fore, we subdivide the considered time interval 7 into a sequence of times to, . . . , tn, tn+1 
,... T, where (e), and (e)„+ı denote the value of a given physical quantity at time tn 
and t„+1, respectively. Suppose that the discrete set of state variables at t, given by 
Man := {Gans SAn, Han Lam) Pant and the local plastic deformation variable co» at 
time t, are known and the time step size At = t„+ı — tn is given. Then, the goal is to 
determine the corresponding fields at time t,,,1 via the algorithmic approximation to the 
weak formulation (3.59) defined as 
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cure bax] f ruvini av = mat 


h,n+l ._ ABSB,n+l — Bn E SB,n A 
G, := 054 M; - fR He, dV + KA? sens j 


Bo 


hynt+1 .— 
Gent. 


ao = 504 Y pes Te RARBG, nti RA Dy nti Varta.) dV A Aa 
IE “(Gh — N) dV + K Qbn+1| s 


p,h 
hn ab X,n+1 — Qb,n a PLA (rh n+l ry) 
GH += ôr? mpn f ii N ng e 
E q SERIES Fi) 
Bo 
(3.69) 
Hereby, a full-discrete definition of the internal dissipation is given by 
h h Sá 
Denes = Vp Ta+l : a + Vf a a (3.70) 


At 


Taking into account small values for the plastic viscosity parameter np, we propose the 
following approximation 


h h h One 0 
j n n 
Tapi” Pt 8 x Jaa fasa) A EAn (3.71) 
such that the internal dissipation can be recast as 
h h ; Qa,n+1 — RAsAmt ZA, 
Dint n+1 := Vp JIa zu a) rea Ne a a+ Vf H? R 7 z (3.72) 


a At 


for practical reasons. 


Note that we apply a staggered scheme for the solution of the multi-field problem, i.e. 
the displacement field along with the plastic and hardening fields {qan+1, Ce. QAn+1 
and Tae sr} the crack phase-field san+ı and the temperature field 04,41 are solved 
successively. For the time integration of the plastic evolution equations, the construction 
of a return mapping algorithm is most crucial. Therefore, we define a trial state as * 


bir = En) Fs (3.73) 


assuming that no further plastic deformation occurs within the time step. Based on this 
trial state, we evaluate the yield criteria (6.37). If ®2 < 0, then the process is purely 


* For the sake of readability, we neglect the labeling of the spatial approximation in the following. 
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elastic and the elastic trial state is the solution. If on the other hand oP. > 0, then the 
trial state is not admissible and a plastic correction is required. Therefore, we apply an 
exponential integration scheme regarding (3.37) which leads to 


OB. 


(Cr exp[ = 2A Fanti Fn] (CU a= > 
Tn+1 


(3.74) 


and 


Asa) = (AS a) exp | = DREN, a]; (3.75) 


respectively. Note that in contrast to standard von Mises plasticity N,+1 Æ Ny and 
||72n+1|| # 1, ie. the plastic correction has to be performed by the Lagrange multiplier 
AP, 1 as well as the components nan+ı which can be obtained by solving the non-linear 
relations 


EN Hp? 
PP — MA =0 and ead — Nani = 0 (3.76) 
a,n+ 


via an internal Newton-Raphson iteration. In addition, the void volume fraction fn+1ı is 
locally calculated by 


fasi = max fo, 1 — (1 — fo)//aetO?}}. (3.77) 


An overview of the return map algorithm is given in Box 1. 
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P E p 
Fri, CR, On41) 1741) Sn(frozen), 0n (frozen) 


Compute trial state: 

bi, = Fna (CR) Fari 

[a Ma] = eig(bi,) 

fn = max{ fo, 1 — (1 — fo)//det [CF], . 

Si, = I] Ao tr 

Tan = HIRVE E (e? — HEN) na @ 


if Jẹ >1 then 

| Pir = 3 (Jna) (SE) — (JE) — 38 (0n — 00)(1 + (JE) 9) gE) 
else 

| Pe = ET (IE — (JE — 38n — 00) (1 + (JE?) 


end 


solve Gy 


3||Taev tr ||? ¿ tr 
| eh, + 2q1 fn cosh Ea - (1+ (mfn)?) =0 


In+ı 
@P — = p 
i. = Otr 7 Trot 
if OF > 0 then 
solve A}, 1, Man+1 


B| Taev n41 |? 3 n D 
el + 2q1 fn cosh [8g] — (1+ (qfn)?) =0 


5 P Po 
Ont — Tny T pAn = 9 


— Na,n+1 = 0 


with Asa)? = (AS tr) exp[ u ZATAR Nan+1] 


o NE S De 
Compute: An Taev n+ (Ag gaa) Pads ac) 


Tn+1 = Tdev,n+1 F Pn+1Jn+1 l 


else 
P 
An+1 E 0 
e e 
Antl Aate 
Tn+1 = Tdev,tr ae PirIn+ı] 
end 


Box 1: Return mapping algorithm. 


3.5 Numerical examples 


In this section we demonstrate the accuracy and performance of the newly developed 
phase-field formulation for thermo-porous ductile fracture. In particular, two examples 
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are considered. Firstly, a benchmark example with experimental data allows us to explore 
the full range of porous thermoplastic material behavior, so that we end up with a com- 
prehensive validated material setting summarized in 3.2. Finally, a second example from 
the third Sandia Fracture Challenge is used to demonstrate the ability of the model to be 
applied to more complex geometries that produce three-dimensional fracture patterns. 


3.5.1 Tensile test 


10 


-10 9 
Figure 3.1: Tensile test. Reference configuration and computational mesh. 


The following benchmark example has been presented in Ambati et al. |7], where experi- 
mental data of steel 1.0553 are used to adjust the ductile fracture model proposed therein. 
For the numerical investigations in this section, we address the following issues: 


I For the present model, physically more comprehensive model, the results have to be 
in a good agreement with the experimental results in terms of hardening, necking, 
crack initialization, propagation and failure as well. 


II For a physical reasonable formulation of ductile fracture, a non-negative dissipation 
of energy must be taken into account. This has to investigated for the elastoplastic 
behavior as well as for the fracture mechanical behavior, cf. (3.40). 


TIT The modeling of the critical fracture energy density plays an important role for the 
ductile fracture behavior. Regarding the definition in (3.22), we have to investigate 
the effect of gep on crack initialization. 
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Elastic parameters 
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Shear modulus u 73255 MPa 
Bulk modulus x 150000 MPa 
Plastic parameters 

Yield stress yo(00) 345.6 MPa 
Ultimate yield stress Yoo (00) 688.25 MPa 
Hardening modulus h(60) 300 MPa 
Saturation exponent wp 16.93 

Thermal softening parameter wh 0.002 K~+ 
Thermal softening parameter wo 0.002 K-! 
Plastic viscosity np 1.1077 MPa-s 
Plastic length scale lp 0.78125 mm 
Initial void fraction fo 0.005 

Gurson fitting parameter qı 1.5 

Gurson fitting parameter qo 1 

Phase-field fracture parameters 

Brittle critical fracture energy Jc., 825 kJ /m? 
Ductile critical fracture energy gep, 175 kJ /m? 
Saturation exponent wf 42.325 
Fracture viscosity ng 1-1077 MPa-s 
Fracture length scale Ir 0.78125 mm 


Thermal parameters 


Specific heat capacity c 


3588 kJ/(m? - K) 


Thermal expansion coefficient 8 1-1075 K7! 
Conductivity Ko 45 W/(m - K) 
Convection Kconv 0W/(m- K) 
Conductivity softening parameter wk 0K-! 
Reference temperature ĝo 293K 
Plastic dissipation factor vp 0.9 

Fracture dissipation factor vr 0.9 


Table 3.2: Material setting of the steel 1.0553 used for the numerical simulations in this 
section. 


IV It is well known for phase-field formulations to brittle fracture, that the fracture 
length scale lp has the character of a material parameter and influences the crack 
initialization significantly. A study of the fracture length scale parameter is indis- 
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pensable, since we suppose that the effect is more pronounced for ductile fracture. 


V To overcome mesh sensitivities we have proposed a gradient extended plasticity 
model as given in (3.17). A comparison of different plastic length scale parameters 
up to a local plasticity model has to be part of a comprehensive investigation of the 
model. 


VI In |7], possible oscillations of the stress field are addressed. With regard to the 
four-field formulation along with the proposed approximation scheme, we have to 
investigate whether oscillations occur within the results. 


VII For the present model, a Gurson-type yield criterion is applied taking into account 
the growth of micro-voids. In comparison to a standard von Mises type yield crite- 
rion, the behavior in terms of hardening, necking and crack initialization has to be 
investigated. 


VIII A temperature-dependent formulation of the hardening behavior is incorporated 
within the present model which determines the elastic response such that we expect 
an impact on crack initialization as well. Therefore, investigations based on isother- 
mal simulations at different temperature levels are useful to understand the model 
in more detail. 


IX For a physically correct representation of the thermomechanical behavior, an energy 
transfer into the thermal field is taken into account via the internal dissipation func- 
tion defined in (3.40). By conducting fully coupled simulations with different values 
of the dissipation factors, the interaction of the different fields can be investigated 
for the entire thermomechanical process. 


X Due to the thermomechanical coupling, the formulation is strongly rate-dependent. 
To investigate this effect, we conduct further simulations and apply different defor- 
mation rates. 


The benchmark problem considers a flat I-shaped specimen of size 140 mm x 20 mm x3 mm 
with reduced width of 12.5mm in the middle. The length of the flaps is 25mm and the 
radius in the transition is 15mm. Figure 3.1 shows the geometry in the reference config- 
uration along with the applied boundary conditions and the computational mesh. The 
outer 20mm of both flaps are subject to Dirichlet boundary conditions. To be specific, 
the lower flap is fixed in space and the upper flap is moved upwards by a displacement 
rate of 0.5 mm/min within a quasi-static simulation setting neglecting inertia effects. For 
the thermal field, no heat in- or outflow is allowed, i.e. the system is adiabatically iso- 
lated. The used computational mesh consists of in total 1936 quadratic NURBS elements 
including a two level local refinement of the region in which plastification, necking and 
crack propagation is expected. 


Regarding point I, Figure 3.2 demonstrates a good agreement of the load deflection result 
obtained by an isothermal simulation at 0 = 293 K and the experimental result given in |7]. 
Therein and in what follows, the displacement is measured at +25mm from the center 
of the specimen. Note that we obtain no significant softening effect due to the crack 
phase-field before crack initialization. In our opinion, this is the correct representation of 
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Figure 3.2: Tensile test. Load deflection result of the present model and the experiment 
given in [7]. 


the problem, since inertia effects are neglected within the simulations. Nevertheless, we 
can reproduce a softer behavior at crack initialization by increasing the fracture viscosity 
parameter nr. 


As shown in Figure 3.3, and in accordance with the experiment, the crack phase-field 
begins to evolve from the center of the specimen, driven by plastic strain localization. 


Next, the non-negative energy dissipation is checked, as addressed in item II. In Fig- 
ure 3.4 we see that both the dissipation due to plastification and the crack growth take 
positive values throughout the process, which is consistent with the second law of thermo- 
dynamics. As for the plastic part, once the yield stress is reached, the system dissipates 
energy, which is enhanced by the strain hardening behavior. In contrast, dissipation into 
the fracture phase-field is much lower until crack propagation. Thereafter, a sudden in- 
crease in fracture mechanical dissipation occurs up to a value of 216.5038 W and plastic 
dissipation ceases as dissipation of elastic energy by crack propagation leads to a lower 
stress level. 


Concerning point ITI, we next consider the effect of different values of the reduced critical 
crack energy density gep. In Figure 3.5, we can observe that crack initiation depends 
strongly on gep. A change of one percent with respect to the initial critical crack energy 
density gee significantly shifts the point at which the material begins to crack, but has no 
effect on the material behavior before that. This is in agreement with the observations 
in Figure 3.2 and our assumption of abrupt failure of steel-like material in a quasi-static 
simulation environment. A similar effect can be observed by varying the fracture length 
scale parameter, as shown in Figure 3.6. Note that this confirms our conjecture in item 


44 3 Phase-Field Modeling of Non-Linear Thermo-Porous Ductile Fracture 


4 
i 1 
| s |- 
E o 
Al 
4 
0.65 
| al 
A r 
A 


Figure 3.3: Tensile test. Result of an isothermal simulation at 0 = 293K. The crack 
phase-field (first row) and the equivalent plastic strain (second row) are plot- 
ted at displacements of u+25 = (10.50, 11.6, 12.5] mm. Additionally, the con- 
stricted cross-section in the middle of specimen is shown, where the solid line 
represents the undeformed configuration. 


IV that the length scale has a large effect on crack initiation and therefore we need to fix 
the value of lẹ for all subsequent simulations. 


Figure 3.7 shows results for different values of the plastic length scale lp, as addressed in 
point V. As expected, the results deviate slightly when necking occurs, but since in this 
case the differences in plastic deformation are too small, there is no visible change with 
respect to crack initiation. Further numerical investigations regarding the plastic length 
scale in |37, 2], have shown the size effects to overcome the unphysical mesh sensitivity 
of the localized plastic deformation and consequently the crack initiation. 


Regarding point VII, an investigation of the modified Gurson model with different values 
of the fitting parameters up to a purely isochoric model with q, = q2 = 0 is shown in 
Figure 3.8 and 3.10. The von Mises stress distribution and the void volume fraction 
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Figure 3.4: Tensile test. Dissipation of energy due to plastification and crack growth. 
The fracture mechanical dissipation part reaches a value of 216.5038 W which 
is not captured in the diagram for reasons of clarity. 
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Figure 3.5: Tensile test. Load deflection results for different values of g.,p. 


at crack initiation are shown in Figure 3.10. Here it can be observed that the necking 
behavior is more pronounced by increasing the Gurson parameters, i.e., by increasing 
the volumetric fraction within the model, which is clearly visible by the deformed cross 
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Figure 3.6: Tensile test. Load deflection results for different values of Ir . 
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Figure 3.7: Tensile test. Load deflection results for different values of lp. 


section in the center of the specimen. This increase in volumetric plastic deformation 
leads to a shift in crack initiation, which is also shown by the load deflection result in 
Figure 3.8. Note that the fracture evolution is governed by the elastic response, so the 
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Figure 3.8: Tensile test. Load deflection result of isothermal simulations at 0 = 293 K 
and different values of the Gurson fitting parameters. 
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Figure 3.9: Tensile test. Load deflection result of isothermal simulations at different 
temperatures. 


maximum stress value at crack initiation is almost identical, but the applied load differs 
due to the different constricted cross sections. 


Concerning point VIII, the temperature dependence of the model is addressed next. 
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Therefore, the load deformation result is depicted in Figure 3.9 for isothermal simulations 
with different temperatures. As expected, higher temperature reduces the yield stress and 
due to this effect on the elastic response, the measured displacement at crack initiation 
is larger, i.e., we obtain larger plastic deformations. This can also be seen in Figure 3.11, 
where the results of equivalent plastic strain and void volume fraction are plotted. 


OvM ovM ovM 
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Figure 3.10: Tensile test. Result of isothermal simulations at 6 = 293K and different 
values of Gurson fitting parameters qı = {0,1.5,1.5} and q2 = {0,1,1.5} 
(from left to right). The von Mises stress (first row) and void volume 
fraction (second row) are plotted at crack initialization. 


Adressing point IX and point X, similar effects can be observed by increasing the 
deformation rate and dissipation factors for fully coupled simulations with an initial tem- 
perature of 293K. As shown in Figure 3.12, an increase in these parameters leads to 
greater and more localized heating, and thus greater plastic deformations, as discussed 
above. However, since the heating is local, unlike the isothermal setting, we obtain a shift 
in crack initiation to smaller displacements measured at +25 mm from the center of the 
specimen, see Figure 3.14. 


Finally, Figure 3.13 shows all relevant physical quantities obtained by the fully coupled 
simulation with ù = 0.5mm/min and vp = vp = 0.9. To be precise, the crack phase- 
field, absolute temperature, equivalent plastic strain, dual strain hardening force, void 
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Figure 3.11: Tensile test. Result of isothermal simulations at temperatures of 6 = {293, 
313, 333} K. The equivalent plastic strain (first row) and the void volume 
fraction (second row) are plotted at crack initialization. 


volume fraction, and von Mises stress distribution at different displacement steps are 
presented. Note that we do not observe any kind of oscillations as addressed in point 
VI. Uniform discretization of all fields, i.e., discretization using quadratic NURBS-based 
shape functions, leads to strong oscillations in the dual hardening force field and stress 
field, causing the simulation to stop. It is evident that the plastic deformation remains 
persistent in case of fracture, so residual stresses remain after the specimen is completely 
fractured. Also, it is evident from Figure 3.13 that the plastic deformation remains 
persistent in case of fracture, so residual stresses are generated within the specimen. 
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Figure 3.12: Tensile test. Absolute temperature result (first row), equivalent plastic 
strain result (second row) and void volume fraction result (third row) of fully 
coupled simulations at crack initialization using different deformation rates 
ù = {0.5, 60,60} mm/min and dissipation factors vr = vp = {0.9, 0.5, 0.9} 
(from left to right). 
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Figure 3.13: Tensile test. Crack phase-field result (first row), absolute tem- 
perature result (second row), equivalent plastic strain result 


(third row), dual hardening force result (fourth row), void vol- 
ume fraction result (fifth row) and von Mises stress result (sixth 
row) at displacements of us =  [7.92, 10.12, 10.94, 11.26, 11.46, 
11.47, 11.51, 11.52, 11.56, 11.58, 11.59, 11.65] mm. 
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Figure 3.14: Tensile test. Load deflection result of fully coupled simulations using 
different values of the deformation rate and dissipation factors. 


3.5.2 Third Sandia Fracture Challenge 


This last example illustrates the applicability of the thermo-porous fracture model to 
complex three-dimensional geometries. Therefore we consider the geometry of the third 
Sandia Fracture Challenget.The geometry in the reference configuration along with the 
applied boundary conditions and the computational mesh is shown in Figure 3.15. The 
examined specimen of size 56 mm x 18 mm x 8mm with a reduced width of 12 mm in the 
middle contains multiple internal channels and cavities. The initial temperature of the 
body is 293 K and no heat in- or outflow is allowed. For an efficient and simple construc- 
tion of the computational mesh based on the given CAD data, the geometry is discretized 
by 111492 Lagrangian tetrahedral elements, whereby a quadratic approximation is used 
for the displacement field, the crack phase-field and the temperature field and a linear 
approximation for both hardening fields. Within the quasi-static simulation setting, the 
lower end of the hollow body is fixed in space while the upper end is moved upwards 
by a displacement rate of 0.2 mm/min Regarding the material setting, steel like parame- 
ters are applied as given in Table 3.2 along with plastic and fracture dissipation factors 
Vp = vr = 0.5 and GTN parameters fo = 0.0005, qı = 0.75 and q2 = 0.5. 


Figure 3.18 shows the crack phase-field, the equivalent plastic strain, the von Mises stress 
distribution and the absolute temperature at different displacements steps. The crack 


t CAD data are taken from the following webpage: 
https: / /drive.google.com/drive/folders/OB8nFwdOwWu0YRkJ4QO0VJWEVCdGM. 
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30 


20] 


Figure 3.15: Third SFC. Reference configuration and computational mesh (sectional 
view). 


phase-field initiates near the intersection between the through hole and the angled chan- 
nels which is driven by the plastic strain evolution. Subsequently, the crack propagates 
along the channels as shown in Figure 3.18 (first row). At the fully broken state, plas- 
tic deformations remain persistent within the specimen inducing pronounced residual 
stresses, see Figure 3.18 (third row). Regarding the thermal field, a global cooling down 
can be observed in the beginning due to the expansion of the specimen. Afterwards, the 
specimen warms up due to plastification and fracture induced energy dissipation. At the 
fully ruptured state, an increase of temperature of about 2K can be observed, see Fig- 
ure 3.18 (fourth row). Note that this warming affects the hardening behavior and thus 
also the fracture behavior as pointed out in the previous example. Moreover, the energy 
dissipation occurs above the conical channel and the internal cavity which impedes an 
heat flux into the lower part of the specimen such that the upper part warms up slightly 
faster. 


Eventually, the load deflection curve is depicted in Figure 3.16. Therein, we can observe a 
pronounced hardening and necking behavior before crack initialization. Figure 3.17 shows 
the results of the crack phase-field at fully ruptured state; the fully broken elements have 
been removed for visualization purpose. 
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Figure 3.16: Third SFC. Load deflection result. 


Figure 3.17: Third SFC. Results of the crack phase-field at fully ruptured state. 
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Figure 3.18: Third SFC. Crack phase-field result (first row), equivalent plastic strain 
result (second row), von Mises stress Bda (third row) and absolute tem- 
perature result (fourth row) at displacements of u = [1.2, 2.3, 3.9, 5] mm. 
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4 Generalized Continuum Mechanics 
for Micro-Structured Materials 


Although most materials look macroscopically homogeneous, on sufficiently small scales 
they usually encompass a distinct microstructure related to their individual microcom- 
ponents, e.g., blood cells, grains, single crystals, fibers, etc. As long as the spatial scales 
of the continuum differ by many order from those of the underlying microstructure, the 
classical continuum approach presented in the previous chapter can account for complex 
material behavior, especially in combination with elaborate constitutive models. How- 
ever, if a structure under consideration is comparatively small, i.e. approaching a typical 
length of the intrinsic microstructure, size effects appear, typically characterized by as 
a stiffer response for a smaller specimen size. In the absence of a characteristic mi- 
crostructural length scale, e.g. size of embedded fibers or spacing between grains or 
crystals, classical continuum theory cannot reproduce these macroscopic manifestations 
of microstructure. 


The main objective of this chapter is to introduce the so-called generalized continuum 
(or microcontinuum) theories to overcome these limitations. Generalized continua can be 
classified into two main groups, see Figure 4.1 and |81, 80, 50]. Higher grade continua 
are characterized by higher order spatial derivatives of the displacement field, such as in 
the second gradient theory of Mindlin [92] or by the gradient of internal variables as in 
strain-gradient plasticity |1, 99]. Higher order continua, generally known as micromorphic 
continua, are represented with additional degrees of freedom that are a priori independent 
from the usual three degrees of freedom of the continuum point [21]. 


After presenting the basic kinematics following |21], of this two main groups of microstruc- 
tured continua, we will focus on the micromorphic theory of Eringen and Mindlin |44, 92] 
for linear elasticity and demonstrate how it incorporates many other models as special 
cases, such as the Cosserat |25] or second gradient theory [78]. The latter will be used 
in the following chapters to account for size effects in the modeling of fiber composite 
reinforcements. 
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4.1 Kinematics 


4.1.1 Configuration and motions 


Higher-order gradient theory 


Bo 


Figure 4.2: Nonlinear deformation mapping 


Considering a three-dimensional body in its Lagrangian and Eulerian configuration Bo 
and B, the deformation of a material point P, identified by its position vector X and a, 
is described via the mapping 

x := (X,t). (4.1) 


We additionally introduce the deformation of a material line L. The material tangent at 
the material finite line element S around the material point P with material placement 
X is denoted by AX, as depicted in Figure 4.2. In the vicinity of P, the following Taylor 
series expansion may be used to generate the mapping of this finite tangent element 


Az := Vy -AX + ZV(Vg) :(AXGAX)+H(AX?). (4.2) 


In the classical Cauchy-Boltzmann continuum, only the linear tangent map Ve is taken 
into account, and thus neglecting the quadratic term V(V«p) and terms of cubic and 
higher degree (6 X+). This requires that the dimensions over which macroscopic gra- 
dients occur are much larger than characteristic length scales of the microstructure. As 
a consequence, only infinitesimal line elements or the tangent vector dX on S at the 
position X affects the deformation map. 


da := Vyp-dx. (4.3) 


When the linear term in the nonlinear deformation mapping (4.2) is combined with the 
quadratic or even higher terms, the classical continuum is extended to a gradient contin- 
uum. When third- and higher-order terms are omitted, this corresponds to a second-order 
gradient theory, which is a specific instance of Mindlin’s second gradient model |91]. As 
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will be explained in more detail later, a size effect is introduced by the second-gradient 
V(V@), which can be related to the curvature in the material line and thus will be used 
to model the size-dependent curvature effects of the embedded fibers. 


Higher gradients of internal fields are used analogously e.g. in theories of gradient plas- 
ticity or gradient electromechanics, see |41, 1, 99, 100, 101]. 


Micromorphic media 


p 


B 
Figure 4.3: Micromorphic deformation maps. 


Rather than considering higher gradients of deformation, the micromorphic theory extends 
the classical continuum in a more general approach by introducing additional degrees of 
freedom per continuum point. More precisely, at each point P in the (macro)continuum, 
a microcontinuum Bo with a set of directors &, with a = 1,..., N is attached, see Figure 
4.3. These directors represent additional degrees of freedom required to characterize the 
microstructural behavior, and denote the orientations and intrinsic deformations of the 
material points of P. Particularly well demonstrated in the physical representations of 
microdeformation in crystalline solids |21], both the geometric point P and the vectors Za 
have their own independent motion. In contrast to highe- order micromorphic continua 
where the directors are higher-order tensors, we focus on first-order micromorphic continua 
with N = 1 and second-order directors. 


In addition to the classical deformation mapping, the so called macromotion 


x := p(X,t) (4.4) 
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we define a corresponding mapping for the micromotion, mapping the attached directors 
from the reference configuration into the current one 


& = p(X, X,t). (4.5) 


In comparison to the macroscopic scales of the body, the material particles are consid- 
ered to be of very small size. Consequently, the micromotion 4.5 is generally linearly 
approximated as 


2:=9(X,t)-X (4.6) 
with the second-order microdeformation tensor p(X, t). 


As originally defined by Eringen |44], a material body is referred to as a micromorphic 
continuum (of grade one), if its motions are described by (4.4) and (4.6) and have unique 
inverses 
X=p9 (Xt 

e (X,t) (4.7) 

X=pU(X,t)-2. 
In order to retain the right-hand screw orientations of the frame of reference, and ensure 
the physical assumption of continuity, indestructibility and impenetrability of matter, we 


assume 
det(Vy) > 0 (4.8) 


and 
det(p) = 1/det(p*) > 0. (4.9) 


Figure 4.4: Deformable directors of the micromorphic continuum 


As shown in Figure 4.4, a material point possesses, in addition to the three usual trans- 
lations, three deformable directors which represent the degrees of freedom arising from 
microdeformations and are denoted by 


(4.10) 
D; = Gig (X,t)E, 
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with the components cartesian unit vectors E, and e;(i = 1,2,3) in the material and 
spatial configuration, respectively. 


4.1.2 Rotation and subclasses 


With the three distinct deformable directors, and its total of 12 degrees of freedom, the 
general micromorphic continuum formulation is appropriate for the modelling of many 
complex materials like polymers with flexible molecules, liquid crystals with side chains, 
animal blood with deformable cells, etc. Its microdeformation can be illustrated with the 
deformable director triad shown in Figure 4.5. 


Moreover, the micromorphic continuum includes many subtheories, each of which is de- 
fined by imposing specific constraints on the microdeformation, see Figure 4.1. A detailed 
systematic presentation can be found in [50]. In what follows, we will briefly discuss the 
most common of these, namely the micropolar, the mircrostretch and the microstrain 
continuum. Later in the chapter, we will show how a second-gradient theory can be 
obtained by imposing analogous kinematic constraints on the generalized micromorphic 


continua. 
R-U 
oe —_—> A 
Bo B 


Figure 4.5: Micromorphic deformation (deformable director triad) 


Just like the deformation gradient introduced in Chapter 2, the microdeformation tensor 
includes both stretch and rotation 


p=R-U, (4.11) 


where R is the orthogonal microrotation tensor (i.e., R-t = RT and det(R) = 1) and U 
is the symmetric and definite positive right microstretch tensor. 


u — R> y 
Bo B 


Figure 4.6: Micropolar deformation (rigid director triad) 
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A micropolar continuum (often simply referred to as Cosserat continuum) as illustrated 
with the rigid director triad in Figure 4.6 may only experience rotation 

#=R-X. (4.12) 


The micropolar theory is typically used for the investigation of granular or multimolecular 
materials, like rigid chopped fibres, elastic solids with rigid granular inclusions, liquid 
crystals, animal blood, etc. |15, 120]. 


Figure 4.7: Microstretch deformation (extensible director triad) 


Besides a rotation R, also one scalar stretch variable A is considered to account for 
isotropic expansion and contraction of a microstretch continuum without microshear- 
ing 


z=AR-X. (4.13) 


In addition to the micropolar director triad with fixed angles, the director triad of is 
isotropically deformable as shown in Figure 4.7. Chopped elastic fibres, animal bones or 
inviscid liquids are typically simulated as a microstrech continuum [50]. 


In a microstrain continuum, the rotation is neglected, rather only the stretch tensor 
U is considered 
£=U-X, (4.14) 


as illustrated in Figure 4.8 . The director triad allows for pure stretch, which is charac- 
terized by the fact that the principle directions maintain their direction. The microstrain 
theory is used e.g. for metallic foams [50]. 


Figure 4.8: Microstrain deformation (stretchable director triad excluding rotation) 
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4.2 Micromorphic linear elasticity 


With the introduced kinematic framework for generalized continuum mechanics, a mature 
and growing field of research has emerged over the last 50 years to conceptualize appro- 
priate constitutive theories and corresponding deformation measures for a wide range 
of microstructured materials. In the following chapters, we will present a specific mi- 
crostructured constitutive formulation at finite strains for the mechanical simulation of 
fiber composite reinforcements. However, since the main objective of this chapter is to 
give a concise introduction to general microcontinuum theories and demonstrate how they 
can be suitably constrained to obtain more specific models, i.e., a second-gradient theory, 
we next present the linear theory of micromorphic elasticity developed by Mindlin [92]. 
Moreover, the paradigmatic structure of the linear model allows a more straightforward 
understanding of the general concept, by avoiding the mathematical complexity brought 
by nonlinearities. This specific linear model has been used to describe a variety of phys- 
ical phenomena, such as polymers with deformable molecules, biological tissues, simple 
harmonic waves, etc. [43]. 


4.2.1 Deformation measures 


Rather than the (macroscopic) deformation map ¢, we consider the classical displacement 


field in the small strain theory 
u=x-X. (4.15) 


Moreover, Mindlin specified the following objective Lagrangian strain measures for in- 
finitesimal deformations, namely the classical linearized macro deformation strain ten- 
sor 


E= (Vu + Vu"), (4.16) 
the relative micro-macro deformation 
y=Vu-%, (4.17) 
and the gradient of the microdeformation (third-order tensor) 


R=VG. (4.18) 


4.2.2 Energetic response 


The general form of the micromorphic strain energy density introduced in |92] is a function 
of the 42 kinematic variables E,;, 7; and Rijk, and is of the type 


1 1 1 
Y = (Ejj, Jij, Rijk) 55 tijk Big En + FLIM Vis Vet + illo Rima 


+ dijkim Vig Rim + EijkimKijk-Eim + fijktij Eri- 


(4.19) 
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Due to the symmetries of the strain tensors, only 903 of the 1764 introduced scalar 
material coefficients are independed. Assuming an isotropic material, the number of 
independent constitutive parameters further decreases to 18 (see |92]), leading to an 
energy function 


1 1 1 1 
v =AEG Es) + Ey Big + ZN + FVV + BVN 
1 1 

+ OR Raji + Ri Raj + Ze iin Risk + g ijj Rik 
1 1 (4.20) 

FeR Rer + JR jihkjR + CT Rijn Rijn + ORAR 

1 1 1 
F Co Risk Ring + 3 C10 Rijn Rin + Rei 
+ fiu Ejj + (Ya + 15) Es. 


Depending on the particular microstructured material behavior under consideration, more 
specific and simplified formulations have been developed. For instance, an energy density 
function with only 6 constitutive coefficients was introduced in |78], which qualitatively 
still reflects all material characteristics of the general micromorphic model of Mindlin 


Y = pellsym(Va — ||? +E (Va — B))? + mallsym(G)|l? + Fra) 
(4.21) 


= a = 
pollskow(Vu - Bl? + IVO. 


Here, the microscopic Lamé moduli u, and A, are related to the response of a represen- 
tative volume element of the substructure. Along with the classical macroscopic Lamé 
moduli uy and A obtained in experiments for large material samples with assumed hetero- 
geneous structure and a characteristic length of zero, the isotropic scale transition param- 
eters u. and A, can be determined from homogenization theory, see [78]. The Cosserat 
couple modulus wu. which governs the asymmetry of the force stresses is i.e. used to 
describe the mechanical behavior of very particular metamaterials as lattice structures 
and phonon crystals. Finally, a, accounts for different microstructural effects, i.e. the 
curvature depencence or some particular optic waves propagation. 


Along with the following restrictions for the introduced material parameters 
Me >0, Pe >O0, 3AC +2ue > 0, pn > 0, 3An+2un > 0, ae > 0. (4.22) 


the well-posedness and positive definiteness of the formulation has been outlined in [78]. A 
thorough identification of the relations for the 6 coefficients proposed here to specific micro 
and macro deformation modes, and the relations in terms of the parameters introduced 
by Mindlin is outlined in [78]. 


Although we will ignore inertia effects in the scope of this work, we briefly introduce the 
micromorphic kinetic energy potential for the sake of completeness, as 


1 Taori 
T= >? lull? + >? lgl? (4.23) 
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with the macroscopic and microscopic mass densities p and p defined per unit of macro 
volume. 


Considering the dimensionless of the second order tensor ¢, the coefficient p has a dimen- 
sion of a bulk density multiplied by the square of the length. With the real density of the 
microstructure po’ , and a characteristic length 1 directly associated to the characteristic 
size of the microscopic inclusions, the microscopic mass density may be written as 


p=, (4.24) 
see |78]. 


Note that in contrast to this special form with a single characteristic length, there are 
also more general types of kinetic energies that take into account different characteristic 
lengths /;; and therefore for more intricate microstructures [92]. 


4.2.3 Strong form 


Ignoring inertia effects, the equations of motion for a linear micromorphic continuum will 
be derived through the principle of virtual work ¿Wi — we = 0. 


For the elastic potential energy given in (4.21), we derive the internal virtual work as 


y y y 
ES :öE + a : 07 + = :-58) dv. (4.25) 


Bo 


Regarding the partial derivatives therein, we introduce the work conjugate stresses to the 
strain components 


Ow 

opm ok 
ow 

als ae (4.26) 
Ow 

S JA 


namely, the Cauchy stress, the relative stress and the double stress, respectively. 


Using the symmetry of a; and the definition of the corresponding deformation measures, 
we obtain 
sw = fa :0Vu+o,:(0Vu-—9p)+6 :-dVeGdV (4.27) 
Bo 


After an integration by parts 


wirt A V-(01+03):0u+V-(du-(0,+07))-(07+V-6):9p+V-(9p : S)dV (4.28) 


Bo 
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the application of the divergence theorem with the unit normal vector N to the boundary 
surface OB, yields 


swt = f Vater) bu (0 +V-6)-5paV+ (0 +0)N)-u+SN :5pdA. 


Bo Bo 
(4.29) 


Note that all further integrations by parts incorporating the boundaries 6? Bo and 8? Bo 
and subsequently leading to external edge and wedge forces are omitted for the sake of 
conciseness. 


Following |78] the external contributions can be formulated as 


wo = [b-öu+ B:öpav + |t-ðudA+ fT: 5704 (4.30) 

Bo T$ PE 
where b is an external body force per unit volume, B an external double force per unit 
volume, t is an external traction force per unit area and T is an external double traction 


force per unit area, acting at the respective Neumann boundaries Tí and TY (see ? for 
further details). 


Eventually, assuming arbitrary variations of the introduced kinematical fields du and 62, 
we obtain the local form of the problem as 


V-(o1+02)+b=0 


(4.31) 
V-G+o,+B=0 
supplemented by the boundary conditions 
u=up on Tŷ 
e or on. Ts (4.32) 


(o, +03N=t on T$ 
RN=T on TẸ 


with prescribed fields up and fp, at the mechanical Dirichlet boundaries I’ and re : 


4.2.4 Special case: second-gradient / strain-gradient theory 


We have already shown in Section 4.1.2, it is possible to impose appropriate constraints 
on the additional degrees of freedom inherent in a micromorphic theory in order to derive 
more specialized models, such as the micropolar/Cosserat theory. In this final section, we 
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will derive the second-gradient model from the micromorphic theory for linear elasticity 
introduced earlier. 


In particular, we have to constrain any relative micro-macro deformation 
y=Vu-p=0 (4.33) 


and subsequently enforcing the microdeformation to be equal to the gradient of (macro) 
displacement 


p+ Vu. (4.34) 
The double stress then reads au 
6 = —_ 4.35 
OVVu’ we) 
leading to the internal virtual work 
owt = fa :6Vu+6:-6VVudV. (4.36) 
Bo 
and the equilibrium condition 
Viloj= 16) Eb=0: (4.37) 


Regarding the simplified strain energy function (4.21), this limit case can be numerically 
obtained e.g. by simultaneously enforcing He > © and ii. —> œ and thus obtaining a 
function of the first- and second-gradient of the displacement 


U(Vu, VVu) = un |sym(Vu)| + tu)? + = [VVul?. (4.38) 


Note that alternative equivalent formulations can be derived with the strain tensor and 
its gradient V((E, VE) or analogously for finite strains with the deformation gradient 
and its gradient W((F, VF) as shown in the next chapter, thereupon "strain-gradient" 
and "second-gradient" theories are often used interchangeably. 


5 A Non-Linear Strain-Gradient 
Formulation for Fiber-Reinforced 
Composites 


In this chapter we focus on a generalized continuum formulation with higher-order contri- 
butions for mechanical simulation of fiber materials. In the modeling of fiber-reinforced 
composites, it is well established to consider the fiber direction in the stored energy in 
order to account for the transverse isotropy of the overall material, induced by a single 
family of fibers. However, this approach does not account for the length scale dependent 
size effects, i.e. the in-plane flexural resistance of the fibers. By using a generalized 
continuum model based on the strain-gradient, the gradient of the fiber direction vector 
can be considered as an additional parameter of the stored energy density function. As a 
result, the improved model presented in this chapter accounts for the bending stiffness of 
the fibers, thereby allowing independent material modeling without recalibration for the 
specific fiber direction or load case. Along with additional material parameters, increased 
continuity requirements on the weak Sobolev space follow in the finite element analysis. 
Here, the isogeometric approach offers a framework that satisfies these weak formulation 
criteria. 

The chapter is structured as follows: Firstly, a second-gradient theory is applied to 
Kirchhoff-Love shell elements. In particular, we embed a model of a woven fabric as pre- 
sented in the work of Steigmann |116] to account for the in- and out-of-plane flexural resis- 
tances of the fibers. Shells, as dimensionally reduced approximations of three-dimensional 
continua, enable a more straightforward presentation including a micromechanical and 
empirical motivation as well as a detailed verification of the higher-gradient contribu- 
tions with analytical results. Moreover, the higher-gradient formulation of the fabric is 
validated with experimental measurements on organic sheets. Finally, a corresponding 
formulation for fiber-reinforced materials will be extended to general three-dimensional 
generalized continua, and verified by the means of two bending tests. 

Note that throughout the chapter plastic, damage, fracture, and thermal effects as well 
as fiber-matrix interactions are omitted for now and considered later. 
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5.1 Shell element formulation 


5.1.1 Micromechanical and experimental motivation 


The classical constitutive equations are homogeneous and do not include any natural 
length scale, and therefore cannot account for any size effects, such as effects of fiber 
diameter or fiber spacing, observed in real fiber composites.To illustrate the theoretical 
basis for expecting size effects, it is sufficient to consider pure bending in plane strain of 
a linear elastic plate with variable Young’s modulus. In terms of rectangular Cartesian 
coordinates (x,y,z), suppose the middle surface of the plate lies in the plane y = 0, 
and the lateral surfaces are y = +h. Assuming that the plate undergoes a pure bending 
deformation, as in the elementary Euler—Bernoulli bending theory, the displacement in 
the x direction is 


vy 
u = —, 5.1 
E (5.1) 
where R is the radius of curvature of the deformed plate. Further, suppose that the 
extension modulus E in the x direction depends on y, so that E = E(y), with mean value 


Eo, where 


h 
1 

Eo =— | E(y)dy. 5.2 
=, di (y)dy (5.2) 

—h 

Then the stress component dy, is 
Ou y 

on = E(y)— = Ey, 5.3 
o. (ya TEOR (5.3) 


and the bending moment applied to a section x = const., per unit length in the z direction, 
is 


h h 
1 
M, = Joey = z J POWA, (5.4) 
—h =h 
and so the bending stiffness B is 
h 
B=M,R= J Ewa (5.5) 
Lh 


On the other hand, if the extension modulus has the constant value Eo the bending 
stiffness is 


h 
2 , 
Be j Eoy?dy = Boh’ (5.6) 
—h 
Clearly, in general B and Bo are not the same. To take a simple example for illustration, 


let 
E(y) = Eo — E, cos = where d = h/N, (5.7) 
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Figure 5.1: DIC measurement of a 25 [mm] specimen in 45° configuration, colors indicate 
local stretches in horizontal direction in [mm/m]. 


and N is an odd integer. Thus d is a measure of the scale of the inhomogeneity of the 
material. Then from (5.5) 
4APhE, d? 6H, 


B=B = Bo + B=, Bı = —— 
o+ 72 o + 173 1= nm 


Bo, (5.8) 


and so the bending stiffness differs from By by a term of order (d/h)?. In the theory of 
classical fiber-reinforced materials, it is explicitly assumed that the fibers are infinitely 
thin and thus infinitely flexible, with the bending stiffness for a single fiber (represented 
by a mathematical curve in the theory) being zero. This clearly corresponds to the limit 
d/h > 0 in (5.7). To relax this assumption while remaining within a continuum theory, it 
is necessary to introduce a length scale into the theory that effectively endows the fibers 
with bending stiffness. 


The in-plane curvature resistance of fibers exhibits the same microstructural size effects. 
The digital image correlation (DIC) of a thermoplastic composite laminate, also known 
as an organic sheet, is shown in Figure 5.1.The additional penalty imposed by the fiber 
in-plane flexure is seen in the deformation patterns of the DIC measurement in compar- 
ison to the straight red line drawn in the figure. The chapter’s proposed computational 
model for advanced mechanical modeling of fiber materials will be validated using the 
aforementioned prototypical composite material. 


5.1.2 Kinematics 


We consider a fiber-reinforced composite as a Kirchhoff-Love shell, as introduced in Sec- 
tion 2.3 with the therein presented configuration and enhanced kinematics. The derived 
standard deformation measures aag and Kag with respect to stretch and out-of-plane 
curvature and the additional measure S%, for in-plane curvature are first examined for 
material symmetry and objectivity. Furthermore, specific deformation measures for the 
embedded fibers are derived and discussed. 


For the matrix material, we use the standard Kirchhoff-Love strain measure, namely the 
right Cauchy-Green tensor and the in-plane Jacobian determinant, as outlined in Section 
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2.3 


and 


Frame invariance 


The present formulation is based on the work of Steigmann |116], who used the Mur- 
doch-Cohen |94] definition of material symmetry to derive the canonical forms of the 
strain-energy functions for woven fabrics with bending resistance. In the following, the 
invariance of the introduced measures {dqg, Ka,3, Sa g} under superposed rigid body mo- 
tions will be briefly demonstrated. This includes translational as well as rotational in- 
variance of the strain energy which in turn implies satisfaction for conservation of linear 
and angular momentum, see Hesch & Betsch [58] and Marsden & Ratiu [79] for further 
details. 


We define rigid motions of the form 


rt :=u+Qr, (5.11) 


where u € E’ and Q € SO(3) is a rotation tensor with the property QTQ = I. Then, 
the first and second partial derivative read 


a =Qra and als = QT og. (5.12) 
From the covariant metric coefficients we obtain 


at, = at: as = (Qr a): (Ar) =T a: QTQT g = aog (5.13) 
=I 
which clearly holds also for the contravariant metric coefficients, i.e. (a°f)" = a. With 
regards to this finding, we can demonstrate 


Stay = Tiaa Tasai = (Qr a) (QT a)i Agaa = 7 QQ rap—Trgdor = Sapa (5-14) 
=I 


as well as 


(Saa) = Skoa la) = Sagra”? = Sip; (5.15) 


Q, 
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i.e. the frame invariance of the co- and contravariant in-plane curvature tensor is con- 
firmed. Eventually, we obtain for the out-of-plane curvature tensor 


Kos = Bas dla 
Qa, x Qaz 
= = Ba Q MIA AA 
Tr Qa x Qa 
_B det(Qr.as, Qaı, Qaz) 
= Bag 
det (az) (5.16) 
det(r a8, @1, 
= Bag — det(Q) Pilas, 21,42) 
= det(aag) 


= Kaf - 


confirming again the frame invariance. 


Fiber deformation measures 


Following Steigmann and dell’Isola [117], we now derive additional deformation measures 
for the embedded fibers which all can be considered as functions of dag, Sig, Fap and 
9%. Thus, we focus on the mathematical modeling of fibers as curves with appropri- 
ate kinematical and constitutive structure to capture the mechanical properties of the 


reinforcement. 


The two fiber threads interlaced at right angles to each other, the so-called warp and weft 
threads of the woven fabric, are characterized by their constant orthogonal vector fields 
in the reference configuration M and L. Described in the reference mid-surface Q, and 
not necessarily with unit length, they can be written as 


L=L°A, and M=M*°A,. (5.17) 


On the deformed mid-surface w they are denoted by l and m and can be represented as 


[=l°a, and m=m°a,.. (5.18) 


Since the fibers are convected as material curves, the tangent vectors are related by 


L M 
àl = F— and Al = F—_ (5.19) 
i ILI í M| 
where we made use of the fiber stretches 
Vaag LLP Vaag LLP Vaag MeMB 
and A ; 5.20 
~= |i- VAT = |ia- am 60 
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The corresponding contravariant components are connected by 
L° = X, J| LI = Ya. LEI" and M® = d3||M||M° = yap M”M”m*. (5.21) 
Furthermore, the change of angle between the fibers can be formulated as 


Aap L“ MP ag Le MP 
P := acos FFT — acos | 22]. (5.22) 
VJ Ap Y Asp MI V apv L” L” aop M° MP 


Hereby, the first term is 7/2 for initially orthogonal fibers. 


To study the change of the tangent vectors l and m within the deformed mid-surface, we 
introduce the auxiliary unit tangent vectors 


p=nxl, q=nxm. (5.23) 


Moreover, let 1: B > E? and m: B > E? be such that 


1(9°) =1(r(0%) and m(0%) = m(r(0°)). (5.24) 
Then the gradients 

Vi(r) =1,@a* and Vm(r)=m,08a" (5.25) 
are easily obtained by inserting d0% = dr - a” into the expressions dl = 1,,dé° = 
laldr - a“) = (la ® a*)dr and dm = m „d0° = (Mma ® a“) dr, respectively. Since 
|| = ||m|| = 1, we have 1-1, = m-m, = 0 and the partial derivatives of the unit 


vectors l and m can be expressed in the right-handed orthonormal bases {l,p,n} and 
{m, q, n} as 
la = (la: p)\p+ (la n)n and m„=(m.:q)lga+(m.:n)n. (5.26) 


Using (5.18) together with the Gauss and Weingarten equation (2.34), the normal com- 
ponents in (5.26) can be expressed using the covariant components of the mid-surface 
curvature bag as 


l.-n= (Pag). n= lbg and Ma'n = (mag) a'n =m’ bga. (5.27) 


The change of the unit tangent vectors in tangential direction can therefore be represented 
by 


Vl(r)l =l la =mp+rn and Vm(r)m = mMm a = Nma + Kmn (5.28) 


with the normal fiber curvatures 


bag Le LP 
yy LE LY 


bag MAME 


„Ze a — Br: 
Ky = dal SS = a MEM? (5.29) 


and Km := bagmém 


and the geodesic curvatures 


m:= la: p and N:=mM"Ma:q. (5.30) 
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The change of the unit tangent vectors of one fiber family in direction of the other fiber 
family leads to 


Vi(r)m = m°l.=dıp+rn and Vm(r)m=m°m, =onqtmn, (5.31) 
where we recognize the torsion and the respective Tchebychev curvatures, |118], 
bag L* MP 
y Qyw L” L” aoo M° MP j 


Similarly to (5.29) and (5.32), we can introduce the normal fiber curvatures kr, Ky and 
the torsion 7 on the reference mid-surface as 


T = bagl"mP = i= m (la: P), bm:= 1% (Ma: q). (5.32) 


Bas L° 1? BasM° MP? i Bas L° M’ 
a A ae Sie peab Tr. oa 
LA A? "MT AMM ONS 9" aA, M M? (5:33) 


Because these definitions rely on unit tangent vectors, the normalization by ||Z|| and 
|M || are involved in here. 


As strain measures that captures the change in the normal curvature of the fibers, i.e. by 
out-of-plane bending, we choose 
Kap LALE 


Kı = A LPL” = Ay — KL and Ka = 


kag MM 
a MAR = Akm — KM - (5.34) 
pu 


The change in torsion by twisting the fibers is given by 
Kop L* MP 
y Ap LELY As Me Me 


Finally, to describe the deformation of the fibers within the tangent planes of the mid- 
surfaces, like in-plane-bending, we introduce as strain measures 


Ks := — = AıAaT OR. (5.35) 


TRIP Scie o MOM? Sate 0 ong Pom L*MP Sapo F 
gL := AnD? > 9m >= A Mim an = ATLA, MONO 


(5.36) 


For initially straight fibers, the deformation measures (5.36) can be represented as 


gL = Aımp+ (L-VAı)l, gu = Aanmq + (M -VA2)m 


5.37 
T= (L $ Vàz)m + A1A20m4 = (M $ VAL + AıAadıp > ) 


with the stretch gradients as VA, = Ay, A® and VA2z = Ag A” (see Equation (58) and 
(59) in |117]). We refer to Section 2.4 of |117], where also the general precurved situation 
is analyzed. However, the restriction to straight fibers lies only in the representation 
(5.37), which should help to get more insight into strain measures (5.36). In fact, these 
strain measures obtained by the contraction of the fiber tangents of the reference mid- 
surface with the third-order tensor a, & S7 represent a combination of the gradients of the 
fiber stretch and in-plane-deformation of the fibers described by geodesic and Tchebychev 
curvatures. 
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5.1.3 Variational formulation 
Energetic response 


For the matrix material we assume the strain energy function per unit reference volume 

Wie (Cy, 6), for the fibers per unit reference area Y. (aag, Sago, Kag, 0%). Moreover, we 

assume that the embedded fibers have a height of h; such that the space occupied by 

matrix material in the reference configuration is defined by {X (0°)|0® € wiso}, where 
h h 


Wiso = [-%, 2] \ (2, t), For the overall composite we obtain for the strain energy 


defined per unit reference area of the surface 


U (aag, Sago, Kap, 0%) = N V= (Ô; (dap, Kap, 9°), 0°) dH’+ Wet (aag, Sago, Kap, 0%) . (5.38) 


Wiso 


Assuming a plane stress distribution in the three-dimensional continuum together with 
a decoupling of in-plane and thickness deformations, we use this specific incompressible 
neo-Hookean strain energy function for the matrix material 


det (Gu) _ 3 


1 5 5 la 
vis — u| tr(C,4G" @ GP +($) -3) =- Q Gor + 
el (Cop PEs PA det (C43) 


0 


) , (5.39) 


with the shear modulus u. Note that this is a correction on a constitutive level for the too 
restrictive kinematical ansatz 2.31, which does not allow for thickness deformation. In the 
numerical examples, we will also consider a compressible neo-Hookean matrix material, 
given by 


g = Fua) 3) 4 aad? 1 — 2In(J)). (5.40) 


Since the plane stress condition cannot be inserted analytically here, we have to condense 
the corresponding conditions numerically. For details we refer to Kiendl et al. [70]. 


Assuming rectangular cross-sections for the embedded fibers, we adopt the following strain 
energy function 


a=1 i=1 

(5.41) 
where the constant material parameters a; are associated to the stretch and the change of 
angle between the fibers. In particular, ag = Eahr correlates to Young’s modulus Ea of the 
particular fiber, whereas a3 = Gh; represents the shear stiffness of a shell continuum. Note 
that the fiber formulation does not take any Poisson effect into account, i.e. the Poisson’s 
ratio y = 0 and thus, G = E in a classical continuum. In the present formulation, 
G represents the stiffness against twist between the fibers and is in general independent 
of the fiber stiffness. The out-of-plane bending stiffness ką depends on the particular 
microstructure of the embedded fibers, i.e. on the second moment of area, such that 
ka = Eala/ba. Here, ba denotes a representative in-plane length scale of the fiber, such 
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that we obtain a bending stiffness of ka = E„h?/12 for a continuously distributed fiber. 
Analogously, we obtain for the in-plane bending stiffness ga = E,b?h;/12. Note that other 
microstructures are possible, although characteristics like fiber height and width have to 
be discussed and validated for the particular geometry. Eventually, the terms associated 
with kz and g3 control the fiber torsion as well as the in-plane deformation captured by 
Tchebychev curvatures and strain-gradients, see (5.37). Assuming that the same fibers 
are used in both directions, the same material parameters can be applied for both fibers, 
i.e. ay = 42, ky = ka and gı = 92. Moreover, we note that the used strain measures of the 
fiber components associated with curvature terms account always for combined effects of 
the gradients of the fiber stretch and bending. Thus, the identification of the material 
parameters in terms of classical beam theory is limited, as we will show in Section 3.5, 


Principle of virtual work 
With the strain energy function presented, we can define the internal energy as 


w™ = l Y dA. (5.42) 
Q 


For the finite element analysis, we discretize the virtual work expression as 
swe = i J BW (Cy; (aag, Kag, 0°), 0°) A? + SW" (das, Sago, Kag, 0%) | AA. (5.43) 


Introducing the stress resultants for the matrix and fiber contributions directly related to 
the mid-surface, respectively as 


Oo gyi» gys» 
nee . , m£ := d me? = 5.44 
Nb ~ Daas > Má Dias and Má Sapo ( ) 
and ayi au: 
nef = / d and mi) := J 20° —— de”, (5.45) 
j Cua í Ôg 


along with Ô; = AON. + 208505) yields 


owt = J @ ma + ne Dias J (mee + mM PY Kap + mee *9Sago | dA. (5.46) 
Q 


Applying the strain energy function of the matrix material given in (5.39) and using 
[Ce#] = [Cag]! together with 


o 


1 1 7 Cro 
x (=>) = -— —det(Ö,)C” = = - —— (5.47) 
0C y \det(Cas) det (Cag)? det (Cag) 
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we get , 
a G = Bra (5.48) 
Cyr 2 det (Cas) 


that has to be integrated in thickness direction to obtain the stress resultants. 
The corresponding stress resultants of the first term (5.44); reads 


na Ov DA, | OV DA) OP ðp awe oco 
fb Oy ðaag Az Odes” 


Op aag OCH daag 
OM Az tan(y) 0p 
A2—1 2 

daag + a2(A2 Band op) daag 

$ (LYLA Syao (LSL Scip) + (MIMAS, 75) (MSM Serp) 
Auer)? OA MEMO 
© (I7 M> So )(LEM" Sey) \ 800% 
OIO ALLA p MOM?) Baap ’ 


=a, (Aı == 1) 
(5.49) 


where we have made use of (5.41) as strain energy function of the fiber material and 


(TATP Bape NI LAS yy) 00 


9L° GL = (A. LPL” > 
(MAME Sago) (M M> Syan) C™ 
. = a 5.50 
JM ` JM (A M" M”) ’ ( ) 
rere (L° MP Sago (L1M*S4y4)C% 
Ay, LY LY As Mo Me : 
Moreover, the relations 
Or 1 Le L£ OA 1 MeM?* 
= P = (5.51) 
daag 2 Vao L L> Apu LELY daag 2 Vao M” MA, M” M” 
and 
Op _ 1 | LM? 
Ga (aag L% MP)? „L# LYa,,M° MP 
f yı - ee NV mae (5.52) 
aag LME 
oo lA LL MO MP + L” L” ao, M? MP 
Wael" Lra, Me MPA Cu we; ) 
as well as the partial derivative of the contravariant metric 
oce oC 
= == C°. 5.53 
daag 0Cag ) 


have been used. For the second term (5.44)2, the stress resultants of the out-of-plane 
bending moments reads 


LeLP Mem? EL? M? 


ech 
ae Aw LED? A MEM JA A, MOM? 


(5.54) 
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whereas the third term (5.44)3 yields 


moh? =g LLP (L L? Soau) gon H ga Me M? (MIM Soyu) oy 
ae (Ap Le L)? (Aj M" M”)? (5.55) 
LMP (LY M> Sayu) fon 
B ALAL A,,M? MP 


The variations of aag, Kag and Sago can be expressed in terms of ôr(0%) = F (0%, eo), 
which is the variational derivative of r* induced by the one-parameter family F(0%, €) for 
which F(9%, co) = r(@%) holds. The variation of the covariant base vectors as well as the 
partial derivatives thereof are da, = Ôr a and aa g = ÔT ag, respectively. The variation 
of the first fundamental form (2.33) can be written as 


lag = Qa bag + Ag ` day . (5.56) 
The variation of the out-of-plane curvature changes is due to (2.51) and (2.35) given by 
Ökap = (90.5: N + aag: ÔN), (5.57) 


where the variation of the unit normal vector reads 


Sn- (I-n®n)(aı x daz + da, X az) (5.58) 


[las x a2|| 


using the identity map I. Eventually, the variation of the in-plane curvature (2.48) gives 
rise to 
OS apo = lab Ay + Qag ` OAs — Dàla -0Qg + 47: 90)). (5.59) 


Strong and weak form 


To identify and collect the resulting bending moments and normal stress contributions, we 
reorganize the virtual work expression. First, we reorganize the variation of the curvature 


tensor as follows 
bbag = N 004.3 + Aug: ÓN , 
dl ea (5.60) 
=n- [aag — Tao] ; 
where we make use of a, - ón = —n- da, and n-ón = 0 along with the Gauss and 
Weingarten relation (2.34);. Taking the second covariant derivative introduced in (2.49) 


into account yields 
dbas =n: [ôT ag Ja S¿300,] : (5.61) 


Next, the covariant term Sago = (cy _ Tà s) gives rise to 
Sapo = apo —Tigdaro , (5.62) 


* For details on the functional space of the admissible test functions see end of Section 5.1.3 
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which yields after some tedious, but straight forward manipulations 
Sago = Ag OT jas + [bagn + Sågar] da, . (5.63) 


Note that we can rewrite the energy equivalently in terms of the contravariant form Sgp, 
which can be rewritten as 


658g = 0” ÖTlap + [bagan — Siga] - day, (5.64) 
see Steigmann |116] for details. The above redefinition of the variation of the different 
strain measures in terms of da, = Ôr ¿ and örjag = ÔT ag — Dogór o allows us now to write 
the internal virtual work as 


Y = N° -ôr „+ M” - drag . (5.65) 


Comparing (5.38) and (5.46) with (5.56), (5.61) and (5.63) yields 


aw \ Yan) IV \ md 
Ne = |2 | — meS sÊ 
Ba) +) Ae 


ow SyM(+A) i ow SYM(yA) ee (5.66) 
ee aj 
and E 5 
ow Ss M(a8) ow S} Maß) 
M** = E i E 
E ) di E ) A en 


Here, the symmetric part of the derivative (e)’””s# is introduced, where the symmetry 
condition is applied with respect to a and ß, ie. 


Ow \ Smer) 1 OV OV 
(2527) (see ® 25.7) : oh) 


3 
Next, we resolve the second covariant derivative 
M! . ôro = M” - ôr ag - MOTS -ÔT a, (5.69) 
such that we obtain for the virtual internal work 
ôv = N° -ôr „+ M” - brag, (5.70) 


where 
N° = N° - MODA. (5.71) 


Insertion in (5.46) and subsequent integration by parts yields 
sw = J (No Me?) ór+ MI- dr] [N -M örda. (6:7) 


„& 


Q 
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The results of the divergence theorem applied on the first term can be decomposed in 
normal and tangential direction, i.e. normal and tangential to the edges of the shell 


‘| [ine - Mi?) -ôr + M” -ôr g| vadS = 
aa 
fer -ÔT )Va — (Mef - ör)va+(M°® - ör „)vavg + (M°® - ör )VaTg dS 


00 
(5.73) 


where v = Va A” is the rightward unit normal to ON and T =7,A” the tangential unit 
vector. Moreover, ôr, = v“ór a and dr, = T*Ór , are the corresponding derivatives, 
decomposed in tangential and normal direction. A secondary integration by parts with 
subsequent application of the divergence theorem of the last term on the right hand side 
of (5.73) on the boundary OQ of the shell yields 


fi M” varg: ôr- dS = Y [M*vargl; - ór; — J (Mv) -ördS , (5.74) 
an en an 


where (e) = d(e)/dS and the square brackets refer to the corner i at the boundary. 
Assuming that the principle of virtual work 0 = 4W** — 5W** is valid with respect to 
the corresponding functional spaces of admissible solution and test functions defined at 
the end of this section, the external contribution can be formulated as 


swe = [g-öraa+ f t-dr+ p-drydS +> fiðri. (5.75) 
Q Y i 


where the edge Y is assumed to be an open set on OQ and the vertices i € = = [1,..,n] 
are defined on 6?Q. Moreover, g denotes a distributed load, e.g. gravitational load. 
Eventually, we obtain the strong form of the second gradient problem as 


(N°-M%),+g=0 on N, (5.76) 


with boundary conditions at the edges 
© ñ (5.77) 


and at the vertices 


(5.78) 


of the shell. As usual for fourth-order boundary value problems, we decompose the 
whole boundary twice. First, Y = Y, U Y;, with respect to Y, N T; = 0, and second, 
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Y = Y, UY, with respect to Y N Y, = Ø. Note that Y rs. are open sets on OQ, and 


we have additionally 2 = =, US, with respect to E; N Z4 = 0. 
Eventually, we introduce the functional space of admissible solutions 
S={renw(Q)|r=FonT,, ry = k on Th, ri =F; on Eg), (5.79) 
and the space of admissible trial or test functions 
V = for € H?(Q) |ör = 0 on Y,, dr, = 0 on Tk, ôr; = 0 on Ey}, (5.80) 


required for the principle of virtual work. Note that the enforcement of the gradient terms 
in the space of admissible solutions is not trivial, see Schuß et al. [108] for details. 


5.1.4 Spatial discretization 


Due to the particular continuity requirements that follow from the incorporation of the 
curvature coefficients into the presented model, the isogeometric method is employed 
for the finite element analysis. Therefore, polynomial approximations of the deformed 
geometry in the shell mid-surface geometry r and its variation ôr are defined as 


rt = 5 R'q, and ör!" = NR q) , (5.81) 


IET JET 


respectively, where q, € R? and dq; € R?. 


The approximations of the tangent vectors and their respective variations read 


a} := r}, = XO Ra and ah := ôr}, = 5 R? ôq - (5.82) 


IET JET 


Thus, the approximations of the normal vector and its variation are given as 


has ai x a a es Me (I -n!® ae: x ba; + dat x al) (5.83) 
laz x az|| ai x az|| 
The approximations of the second covariant derivative and its variation read 
an. g — Tg = > RaQ and dans = ôr? g = No RY 369, - (5.84) 


TET JET 


Moreover, the membrane strains, the in- and out-of-plane curvature as well as their vari- 
ations are discretized as 


ahy =ah-a and daiz,=ah-da}+a}- dah, (5.85) 


Ku = aus -n> and rig = Jans nt ans - in", (5.86) 
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Sao = ah g aÈ Pa, and Sig, = da, sa, + ahg dag == a (5.87) 
To complete the approximation of the strain measures, the discrete right Cauchy-Green 
tensor and its variation are given as follows 

} 3,h h Ah h 35,1 
C= =a ag ~ 20 ang n and Cag = faas + 20%0k, y (5.88) 


Now, using above approximations the discrete version of the internal virtual work given 
in (5.46) reads 


éwineh = I [ing + Tap, a + (m; iso, on + Ms el B T mann OS dA, (5.89) 


Qh 


where discrete versions of the stress resultants for the matrix material are given as 


Ovyisoh Ch Oviso,h feu 
= | (Co) ag? and mea, =f 26° CH) age (5.90) 
zay ach 


and for the fiber material as 
aß UND (an Rs Kig) 


sph = AE j 
Oars 
fib,h h Kh 
apo _— ow (aa a, aßo? Ka) 
Mib,h , (5.91) 
T 
fibh/ h ch „h 
meee _ Ow (Caen Sabo» Kap) 
fib,h ~~ h è 
OS Bo 
Moreover, the discrete version of the external virtual work given in (5.75) reads 
goth — J g- ôr”dA + J t- ôr” + p- (verre, dS +Y fi- dr? , (5.92) 
Qh yh i 
where v? = v®hA, is the discrete version of the rightward unit normal to 00°. 
5.1.5 Numerical examples 
0.1 ] 
Incompressible matrix material 
Shear modulus p 2 [MPa] 
0.05 4 Fiber material 
Young's modulus E 6 [MPa] 
Thickness hr 0.001 [m] 
e | AA AAA A 
fm] 0 0.05 01 0.15 Table 5.1: Material setting of the composite material. 


Figure 5.2: Computational mesh. 
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The objective of this section is to present a series of numerical examples and to demon- 
strate the accuracy and performance of the proposed second gradient model. To begin, 
we examine a series of benchmark examples in terms of analytical solutions, verifying the 
comprehensive fiber material formulation. 


ul 
pr 


Figure 5.3: Tensile test. Problem setting. The lines illustrate the fiber structure. 


Secondly, the combined matrix and fiber material is calibrated and validated with various 
experimental tests under different loading conditions and fiber orientations. 


| =: 
3.8 [kJ /m?] 5.6 


Figure 5.4: Tensile test. Strain energy density of the matrix material at displacement 
u = 0.15 [m]. 


With the adjusted set of material parameters obtained in the previous examples, we ex- 
amine a final example, demonstrating the application of the model to a more complex, 
non-planar geometry. In addition, the influence of the novel in-plane bending contribu- 
tions of the fiber materials is particularly emphasized. 


Verification 


In a step by step verification of the numerical framework, we conduct several numerical 
tests which allow us to investigate each effect of the fiber material separately. Unless 
otherwise defined, a shell of size 1x b = 0.15 [m] x 0.1 [m] discretized by 3 x 2 cubic 
B-spline based elements is applied. The computational mesh and the material setting are 
given in Figure 5.2 and Table 5.1, respectively. Note that for each B-spline based mesh 
considered in this section, knot vectors of the structure [0f =... = 0% 11 <...< Op 41 = 
= 0¢ .,,41] are applied. 

Tensile test 

We start our investigation with a tensile test to verify the behavior related to a stretching 
of the fibers, i.e. we set a, = ag = heE = 12 [kN/m] and neglect all other contributions 
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Figure 5.5: Tensile test. Strain energy density of the fiber material (left) and the 
composite material (right) at displacement u = 0.15 [m] and fiber configura- 
tions of a = [0°, 15°, 30°, 45°] (from top to bottom. The lines illustrate the 
deformed fiber structure. 


of the fiber material. The left edge of a shell is fixed and the right edge is moved in 
e¡-direction by u, see Figure 5.3 for details. 


In Figure 5.4, the strain energy per unit reference area is depicted for the pure matrix 
material using a shell thickness of h = 0.002 [m]. Moreover, strain energy densities for 
the pure fiber material with hp = 0.002 [m] as well as for the composite material using 
the setting from Table 5.1 are shown in Figure 5.5. Therein fiber orientations of a = 
[0°, 15°, 30°, 45°] are applied. Note that the fiber material in the 0° configuration does not 
experience lateral strains under longitudinal displacement conditions at the right edge. 
Here, we obtain a constant energy density of 6 [kJ/m?] which is in accordance with the 
analytical solution, i.e. X> aa(àa — 1)?/2 with Ay = 2 and Ag = 1. 

Figure 5.6 shows the load deflection result, where results for the pure matrix material, 
the pure fiber material and the composite material are considered. Concerning the fiber 
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Figure 5.6: Tensile test. Load deflection result of the fiber and matrix material (left) 
and the composite material (right). 
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Figure 5.7: Shear test. Problem setting. The lines illustrate the fiber structure. 


material with a 0° fiber orientation, we obtain a load of 1.2 [kN] which matches again with 
the analytical solution. This can be calculated in a straightforward manner using (5.66) 
and (5.77) by which we obtain a constant traction at the right boundary of t = 4*e;. 


Shear test 


Next, we deal with a simple shear test. To be specific, the lower edge of a shell is fixed, the 
left and right edges are fixed in es-direction and the upper edge is moved in e,-direction 
by u as illustrated in Figure 5.7. To verify the implementation, numerical results using 
pure fiber material (hr = 0.002 [m]) have to match analytical solutions. Therefore, we set 


a = 0°, aı = ag = 12[kN] and az = htu = 4[kN/m]. All other parameters are set equal 
to zero. 


Figure 5.8 and Figure 5.9 show the distribution of the strain energy density and the load 
deflection result, respectively. Therein, we obtain a homogeneous strain energy density of 
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Figure 5.8: Shear test. 
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Figure 5.9: Shear test. Load deflection result. Numerical and analytical solutions of 


resultant forces at the right edge Yı (left) and upper edge Yə (right) are 
depicted. 


Figure 5.10: Out-of-plane bending test. Problem setting. The lines illustrate the 
fiber structure. 
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3029.3 [J/m?]. The obtained results are in accordance to analytical solutions, which will 
be derived next. 


Assuming a linear parametrization of the reference mid-surface 


r=0'e,+0e, with 9'e[0,1] and & €(0,)] (5.93) 


the analytical solution for pure fiber material with a fiber configuration of a = 0° can 
be derived as follows. Note that we obtain A; = e, and Aa = ez due to (5.93). The 
isochoric deformed configuration of the mid-surface shown in Figure 5.8 is described by 


1 2 1 Pu 2 
r(@,0°) = (0 ars eı +0, (5.94) 


where u is the predefined displacement of the upper boundary, see Figure 5.7. Now, 
applying (2.32) and (2.33) we obtain 


ole 


2 
M=1, A= (5) +1 and p =~ — acos (5.95) 


b 2 uy? 
(=) +1 
b 


as fiber stretches and the change of angle between the fibers, respectively. Insertion of 
the strain measures into (5.41) and (5.77) yields 


2 


2 
2 u 
| (=) +1-1) + Ztan Z — acos Je (5.96) 
| 
b 


and after some straightforward calculations 


0 5 1 
1 II 
N! =al |1 and N?=a | 1-— > LT ul, (5.97) 
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tan ie acos b 
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is applied. Eventually, using (5.77) we can evaluate the resultant forces at the right 
boundary Y, and upper boundary Y, as 


Fy, = f [m n?) vr dS = Nb and Fr, = f [m n?) vr, dS = N?L, (5.99) 


Yı Ya 


respectively, where vy, = [1,0]? and vy, = [0, 1]* are corresponding unit normal vectors. 
Thus, for the final deformation state, i.e. u/b = 1, we obtain 


Yf» (u = 0.1 [m]) ~ 3029.3 3 (5.100) 


as homogeneous strain energy distribution and the forces at the boundaries are 


| 0 | ee 
Fr, (u=0.1[m)) = [800 | [N] and Fr,(u=0.1(m)) = |1127.2| [N]. (5.101) 
0 0 


Out-of-plane bending test 
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Figure 5.11: Out-of-plane bending test. Configuration of the shell and strain energy 
distribution at different load steps. Pure matrix material (left) and pure 
fiber material (right) are applied. The lines illustrate the deformed fiber 
structure. 


For the out-of-plane bending test, the left edge of a shell is clamped and the right edge 
is subject to an external out-of-plane torque u = Mes as shown in Figure 5.10. In 
contrast to previous tests, the shell used for this test has a length of l = 1[m] and a 
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Figure 5.12: Out-of-plane bending test. Convergence results. 


width of b = 0.1 [m] in the reference configuration, see Figure 5.11. Using a thickness of 
h = 0.002 [m] for pure matrix material and a thickness of hp = 0.002 [m] for pure fiber 
material, the area moment of inertia is obtained by I = bh? /12 = bh} /12 = 66.6667 [mm‘]. 
The applied torque M = 27 EI /l is chosen such that the shell describes a perfect circle in 
the asymptotic limit (see e.g. |69]) which holds for the fiber material by setting a = 0°, 
a, = hE = 12[kN/m] and kı = EI/b = 0.004 [Nm]. Moreover, for the matrix material 
we apply the compressible model given in (5.40) with u = E/2 = 3 [MPa] and k = E/3 = 
2 [MPa] which corresponds to a Poisson ratio of y = 0. 
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Figure 5.13: In-plane bending test. Problem setting. The lines illustrate the fiber 
structure. 


In Figure 5.11, the configuration of the shell is plotted along with the strain energy 
distribution for pure matrix and pure fiber material at different load steps. The results 
shown therein are obtained by 80 x 8 cubic B-spline based elements. Moreover, Figure 
5.12 shows the convergence properties for both materials, where the error is determined 
as gap between the ends of the shell. Here, we observe a limit of convergence at a value 
of approximately 1074 [m] for pure matrix material, whereas for pure fiber material an 
error of less than 2 - 1078 [m] is obtained. Note that the limit of convergence using pure 
matrix material can be justified by known locking effects of Kirchhoff-Love shell elements 
as it may occur even for high polynomial approximations, see e.g. |66, 12, 54, 42, 14]. In 
addition, the usage of an internal Newton-Raphson iteration within the implementation of 
a compressible neo-Hookean model for the matrix material also limits the convergence. 


In-plane bending test 
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Figure 5.14: In-plane bending test. Configuration of the shell and strain energy distri- 
bution at different load steps using pure fiber material. The lines illustrate 


the deformed fiber structure. 
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Figure 5.15: Torsion test. Problem setting. The lines illustrate the fiber structure. 


Next, a verification of the in-plane bending behavior of the fiber material is demonstrated. 
Therefore, we use the shell mesh defined for the out-of-plane bending test and set a = 0°, 
a, = hE = 12[kN/m] with hs = 0.002 [m] and gı = 9 = El/b = 10 [Nm] with I = 
heb? /12 = 166667 [mm*]. All other parameters of the fiber material are set equal to zero. 
Moreover, we apply u = Mez as an external in-plane bending torque, see Figure 5.13 for 
details. Again, the applied torque M = 27 El /l is chosen such that the shell describes a 
perfect circle in the asymptotic limit, which is demonstrated in Figure 5.14. Even if this 
result is not physically reasonable for a real composite due to the material penetration, 
it verifies the implementation related to in-plane bending of pure fiber material. 


Torsion test 
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Figure 5.16: Torsion test. Strain energy distribution of the fiber material at Ad = 180° 
using different fiber configurations a = [0°, 15°, 30°, 45°] (from left to right 
and top to bottom). The lines illustrate the deformed fiber structure. 


Eventually, we verify the numerical framework related to a torsional deformation of the 
fiber material, see Figure 5.15. Therefore, the deformation is predetermined such that the 
shell undergoes a constant twist of 96/OX} = 1200 [°/m] with X, = X -eı. Moreover, the 
parameters of the fiber material are defined by hr = 0.002 [m] and kz = uI,/b = 5.002 [Nm] 
with a polar moment of inertia I, = bhg(b? + h?)/12 = 166733 [mm*]. Again, all other 
parameters are set equal to zero. 


Strain energy densities obtained at Ab = 180° are depicted in Figure 5.16 for different 
fiber configurations. As expected, we obtain a highly stressed region along the central axis 
of the shell especially in the 0° fiber configuration, whereas the fiber material does not 
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[/m] 


Figure 5.17: Torsion test. Twist of fiber at Ad = 180° using a fiber configuration 
of a = 0° (left) and 45° (right). The lines illustrate the deformed fiber 
structure. 


contribute to the strain energy in the 45° configuration since the fibers are not twisted. 
This is also evident in Figure 5.17, where the twist of fiber is plotted for fiber configurations 
of a = 0° and a = 45°. Concerning the 0° configuration, fibers in e¡-direction undergo a 
twist of 09/0X, = 1200 [°/m] at the central axis of the shell. Moreover, analytical details 
related to the torsional test will be detailed next, verifying the numerical results shown 
in Figure 5.16. 


Assuming a linear parametrization of the reference mid-surface, cf. (5.93), with 6! € [0, 1] 
and 6? € [—b/2,b/2] such that A, = eı and Az = e», strain energy distributions for 
pure fiber material with different fiber configurations can be evaluated as follows. The 
predefined deformation of the mid-surface shown in Figure 5.16 and 5.17, respectively, is 
described by 


9! 9! 
r(0',0°) = Hleı + 0° cos ES es + 0° sin ES ez. (5.102) 


Applying (2.32), (2.35)2 and (2.51)2 we obtain 


T 
kil = k2 = 0 and kiz =k = Ga; i (5.103) 


Concerning different fiber configurations defined by a, the vector fields L and M in the 
reference configuration are given by 


L' = cos(a), L? =sin(a), M' = —sin(a), M?’ = cos(a). (5.104) 
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Eventually, insertion of these derivations into (5.41) yields 


2 


(Ka = 8 a (cos(a)? — sin(a)?)?. (5.105) 


yib 
2 (Om? +2 


Note that for the evaluation all material parameters are set equal to zero excepted k3. 
Obviously, for the a = 45° fiber configuration, the strain energy distribution related to a 
torsional deformation of the shell is identical to zero, whereas for each other configuration, 
the highest value is obtained at the central axis of the shell, i.e. at 6? = 0. We obtain for 
the a = 0° fiber configuration 


J J 
Yë» (6? = 0) = 1097.06 -z and Yi» (92 = b/2) = Ye (8? = —b/2) ~ 523.25 a 
(5.106) 


for the a = 15° fiber configuration 
J J 
yf» (92 = 0) = 822.8 5 and y (9? = b/2) = YL (0? = —b/2) = 392.44 = (5.107) 


and for the a = 30° fiber configuration 


J J 
(P? = 0) & 27427, and vb (9? = b/2) = V» (8? = —b/2) ~ 130.81 = 
(5.108) 


Calibration and validation 


We investigate Tepex®dynalite 102-RG600(1)/47 from LANXESS as prototypical com- 
posite material, consisting of 47% vol. of a woven fabric (roving glass) and polycapro- 
lactam (PA 6) as matrix material. The matrix material is hydrophilic with a maximum 
absorption of moisture between 2.6% and 3.4%. For the investigations presented here, a 
single layer material with layer thickness of 0.5 [mm] has been specifically manufactured 
for this test, since multi-layered material with different fiber angles intermix various phys- 
ical effects within the measurements. 


Figure 5.18: Tension test. Clamping devise with 36 mm specimen, undeformed config- 
uration. 


The surrounding matrix material stiffens the woven fabric and connects the junctions of 
the fibers. The fibers within the composite material have significant bending stiffness 
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due to the matrix material surrounding every single string within the fibers. Additional 
stiffening mechanisms are present depending on the actual weave of the woven fabric, since 
the fibers are subject to frictional behavior at the junctions. These effects arise only for 
the combined composite material and thus, cannot be separated within the experimental 
investigations, i.e. we always obtain a mean value of the material characteristics. 
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Figure 5.19: Tension test. Experimental force-displacement curves for 25 [mm] speci- 
men with 30° fiber orientation. 


The matrix material is in general governed by a far more complex constitutive behavior in 
the inelastic regime. Moreover, fibers and matrix might separate in the presence of large 
deformations, i.e. both materials can be characterized by separately assigned deformation 
fields, coupled via corresponding constitutive laws. For now, we restrict ourself to the 
hyperelastic regime without fiber separation. 


Tension test 


We first elaborated tension tests with 0°, 15°, 30° and 45° fiber orientation. Here, experi- 
mental data of 25 [mm] specimen are compared with numerical results using 34 x 10 cubic 
elements. Shear and bulk modulus used within the compressible model (5.40) for the PA 
6 matrix material are assumed to be u = 384.62 |N/mm?] and x = 833.33 |N/mm?|, 
respectively. This setting corresponds to Young’s modulus of Ej, = 1000 [N/mm?] and 
Poisson ratio of y = 0.3. Moreover, for the glass fibers we assume Eg, = 73000 [N/mm?]. 
The thickness of the material is h = 0.5 [mm], with 47% vol. of woven fabric, from which 
50% are aligned in L and M direction, respectively. Thus, we obtain the material con- 
stants for the fiber stretch contributions via 0.5hEg, x 0.47 as aıya = 8577.5 [N/mm], 
which corresponds to a combined Young’s modulus of 17685 [N/mm?]. As can be seen in 
Figure 5.20 for the 0° fiber orientation, the simulation fits perfectly the experiments. 


The value for a3 = 250 [N/mm] is related to the interaction between the fibers and can 
not be derived from other constants. Thus, we have fitted the value to the experimental 
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Figure 5.20: Tension test. Force-displacement curves of the tension test for 25 [mm] 
specimen with 0°, 15°, 30° and 45° fiber orientation (left to right, top to 
down). Simulation results are compared with experimental investigations. 
Additionally, simulation results without in-plane bending stiffness (Sim. 
wb) are shown to highlight the global effect of the non-classical in-plane 
curvature term. 


observation obtained by a fiber orientation of 45°, see Figure 5.20 (right, bottom). The 
bending stiffness gı = ga on the other hand is related to Young’s modulus via gij2 = 
Esp I /b = 893500 [Nmm], where the second moment of inertia is / = b’hr/12 with hr = 
0.47h. Note that this approach is only a motivation for choosing the material constants 
gi and ga, since this allows us a comparison with analytical results as shown in Section 
5.1.5. The values k1/2/3 do not contribute in these tension tests and g3 was set to zero. 
Additional simulation results presented in Figure 5.20 demonstrate the effects of the in- 
plane bending term, i.e. we set gı = 92 = 0, to compare the results with an anisotropic 
Kirchhoff-Love shell formulation neglecting the higher-gradient terms. 


This example shows that the bending stiffness has considerable effects on the global 
stiffness and is not negligible. The constitutive parameters chosen with respect to the 
different matrix and fiber materials and the respective geometry allow for a nearly perfect 
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Figure 5.21: Bending test. Deformed configuration. 
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Figure 5.22: Bending test. Force-displacement curves of the out-of-plane bending test 
for a 25 [mm] specimen with 0° fiber orientation. Simulation results are 
compared with experimental investigations. 


prediction of the macroscopic behavior of this material with microstructure. Only the 
parameter az is fitted, since the shear stiffness of the fibers depends on the weave of the 
fabric and the production process in terms of applied pressure, heat and cooling rate. 
This result illustrates the necessity to include higher-gradient terms in the model for 


98 5 A Non-Linear Strain-Gradient Formulation for Fiber-Reinforced Composites 


Figure 5.23: Bending test. Permanent deformation of the sample after a full loading 
cycle. 


fiber-reinforced polymers. 
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Figure 5.24: Tension test. DIC measurement of a 25 [mm] specimen in 45° configura- 
tion (top) and simulation results (bottom) at 4.6 [mm] total displacement, 
colors indicate local stretches in horizontal direction in [mm/m]. Black lines 
indicate the fiber direction in the current configuration. 
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Bending test 


Next, we investigate the out-of-plane bending stiffness, hence we validate the values for 
kı and ky. Therefore, a 4-point bending test as shown in Figure 5.21 is conducted. For 
the experiments, a stack of 5 specimen with 125 [mm] x 25 [mm] in 0° fiber configuration 
were used, such that we expect a stiffness of 5 times kı;2 = Egpl/b = 78.949 [Nmn], 
where / = bh3/12 and hę = 0.47h. The support structure has a width of 40 [mm] between 
the inner brackets, and 30 [mm] between the outer and the inner brackets, such that the 
whole device has a width of 100 [mm]. 


For the numerical validation, we discretize the shell with 50 x 10 cubic elements with 
prescribed displacements in upward direction at the contact points of the support struc- 
ture, such that the shell can slide over the bracket points and the system is subject to 
pure bending. The results in Figure 5.22 demonstrate the consistency of the formulation 
with the experimental observations until a displacement of about 15 [mm]. Afterwards, 
inelastic effects can be observed, resulting in a permanent deformation of the sample as 
can be observed in Figure 5.23. After a displacement of 25 [mm], the sample starts to 
slip over the brackets. 


Figure 5.25: Tension test. DIC measurement of a 36 [mm] specimen in 60° configura- 
tion (top) and simulation results (bottom) at 4.6 [mm] total displacement, 
colors indicate local stretches in horizontal direction in [mm/m]. 


Digital Image correlation 


Finally, we show in Figure 5.24 (top) the results of the Digital image correlation for the 25 
[mm] specimen with 45° fiber orientation and a total displacement on the left boundary 
of 4.6 [mm]. Note that both sides are clamped in such a way that the displacement is 
prescribed and, additionally, the higher-order boundary conditions are also set, preventing 
in-plane twist and torsion of the fibers. A direct comparison with simulation results is not 
feasible, since the current model does not incorporate elastic softening or inelastic effects of 
the polymeric matrix material. However, we want to highlight the strain-gradient effects 
in terms of the resulting deformation patterns. Therefore, the stiffness of the matrix 
material Fis, and the bending stiffness gı and ga is reduced by a factor of 500 to fit the 
globally reduced stiffness of the inelastic deformed composite. The bending stiffness of 
the pure fiber is nearly zero, but contribute considerably to the global stiffness if the fibers 
are embedded within the matrix material. If the matrix is subject to inelastic effects, the 
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bending stiffness is effected as well. Note that we have observed pronounced inelastic 
deformations by evaluating hysteresis curves. A detailed investigation of inelastic effects 
including plastification of the matrix material will follow in the next chapter. In both 
plots in Figure 5.24, the strain-gradient effects result in the stiffness of the fibers against 
in-plane bending, which can be observed at the crosspoint of the fibers as highlighted by 
the red line. 


Using the same material setting, we compare next results for a 36 [mm] specimen in 
30/60° configuration, see Figure 5.25. As can be seen, the deformation patterns fit in 
general, although the fibers at the vertices of the sample already start to fracture. 


As final result of this section, we summarize the determined and validated material setting 
of our prototypical composite material Tepex®dynalite 102-RG600(1)/47 in Table 5.2. 


Warp and weft direction are given in 0° and 90° direction, respectively. Both exhibits a 
slightly different material behavior and are not perfectly aligned due to the production 
process. This can be calibrated using different values for ay/2 and g1/2, which we omit 
here since the general behavior of this second-gradient material is the key issue in this 


paper. 


Compressible matrix material (PA 6) 


Shear modulus p 384.62 [N/mm?] 
Bulk modulus x 833.33 [N/mm?] 
Fiber material (roving glass) 

Tensile stiffness a1, as 8577.5 [N/mm] 
Shear stiffness a3 250 [N/mm] 
In-plane bending stiffness g1, go 893500 [Nmm] 
Out-of-plane bending stiffness kı, ka 78.949 [Nmm] 


Table 5.2: Material setting of Tepex®dynalite 102-RG600(1) /47. 


Shaft example 


This last example demonstrates the applicability of the proposed model to complex ge- 
ometries with non-planar reference configurations. Therefore we consider a shaft of length 
50 [mm] and radius 50/7 [mm] as shown in Figure 5.26, using the validated material set- 
ting given in Table 5.2. Additionally, we set gą = 0 and kz = 0. The discretization of 
the shaft consists of 200 cubic B-spline based elements. To obtain a closed cylinder, the 
extended Mortar method |108] is applied between the opposing edges in peripheral direc- 
tion of the shell, maintaining G*-contimuity, see also Dittmann et al. |40] for more details. 
The displacement on the left in es-direction of the cylinder is fixed in space, whereas the 
displacement on the right is fixed on the (eı, ez)-plane and rotated as illustrated in Figure 
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5.26. Note that no higher-order boundary conditions according to (5.77)3 and (5.78) are 
applied. 
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Figure 5.26: Shaft example. Problem setting and computational mesh. 


First, we investigate a pure matrix material. After a prescribed angle of Ad = 3.3", the 
shell starts to buckle. The last state of equilibrium is obtained at an angle of Ad = 
12.6°. Note that e.g. an arc length method can be used to calculate further states in the 
geometrically unstable regime. In Figure 5.28, the applied torque necessary to obtain the 
prescribed twist at the right hand side of the shell versus the angle of twist Ad is plotted 
for the pure matrix as well as for the composite material with different fiber orientations 
with and without the in-plane bending terms. 
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Figure 5.27: Shaft example. Strain energy density of the matrix material. 


Next, the composite material is investigated. In Figure 5.30, top left, the strain energy for 
a 0° fiber configuration after a prescribed twist of Ad = 90° is displayed. Note that the 
fiber orientation and not the computational mesh is shown. The corresponding torque is 
given again in Figure 5.28. Several issues can be highlighted here. First, the higher-order 
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boundary conditions as presented in (5.77)3 are not set and thus, the fibers can twist in- 
plane at the boundary. This can be observed very well in Figure 5.30, top left. If bending 
would be prohibited at the boundary, the fibers have to remain orthogonal with regard to 
the boundary. Since these additional boundary conditions are unknown to a first gradient 
theory or even to the classical Kirchhoff-Love shell theory which only considers the out- 
and not the in-plane contributions, the proposed formulation is the only way to enforce 
this consistently within the mechanical formulation. Second, the in-plane terms stabilize 
the shell against buckling. Figure 5.31, top left, shows the results if the additional in- 
plane terms are not present, i.e. the material constants gı and ga are set to zero. Thus, 
we obtain the results for an anisotropic Kirchhoff-Love formulation, usual fitted via a 
first-gradient homogenization process, showing large buckling modes. To illustrate this 
further, the strain energy of the in-plane contributions for the deformed configuration (cf. 
Figure 5.31, top left) is plotted in Figure 5.29. As can be observed, a dramatic increase 
of the in-plane strain energy occurs for this heavily deformed configuration, hence, in the 
mechanical equilibrium this buckling mode is not possible. 


Additionally, results for 15°, 30° and 45° fiber orientation are plotted (top right, lower 
left and lower right picture in Figure 5.30 and 5.31). The corresponding torque is plotted 
again in Figure 5.28. Notably, a nearly constant strain energy distribution can be observed 
in Figure 5.30, lower right and top left. 
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Figure 5.28: Shaft example. Torque versus angle of twist. Displayed are the results 
for the pure matrix material, the composite and the composite without the 
additional in-plane contributions (composite*). 


5.1 Shell element formulation 103 


| U A 
0 [kJ/m?] 750 


Figure 5.29: Shaft example. Estimated strain energy density of the in-plane contribu- 
tion using the deformed configuration in 5.31, top left. 
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Figure 5.30: Shaft example. Strain energy density of the composite material with 
different fiber orientations a = [0°, 15°, 30°, 45°] (from left to right and top 
to bottom) plotted at the final deformation state marked in Figure 5.28 
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Figure 5.31: Shaft example. Strain energy density of the composite material with 
different fiber orientations a = [0°, 15°, 30°, 45°] (from left to right and top 
to bottom) without in-plane contributions plotted at the final deformation 
state marked in Figure 5.28 


5.2 Solid element formulation 


In the previous section, we used experimental results to demonstrate that our shell for- 
mulation with higher-order in-plane bending contributions, in contrast to the classical 
Cauchy continuum, allows for independent modeling of the composite without the need 
to recalibrate the material setting for specific fiber orientation. In this section, this gen- 
eralized continuum formulation for long fiber materials is adapted to a more general 
three-dimensional solid element, in a more straightforward manner. As in the previous 
section, we focus on modeling of the mechanical behavior of the fibers, and use a stan- 
dard elastic matrix material. To avoid being too repetitive, we will skip over the spatial 
discretization, which will be discussed in depth in the following chapter. The implemen- 
tation is then verified by numerical examples, focusing on the validating of the bending 
contributions. 
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5.2.1 Configuration and kinematics 


As introduced in Chapter 2 in the context of classical continuum mechanics, we consider a 
body of interest By € R? in its reference configuration along with the deformation mapping 
p(X,t). The deformation gradient F(X,t) with J:=det[F] > 0 can also be defined in 
terms of the principal stretches A, with a = {1,2,3} and the principal directions of the 
left and right stretch tensors n, and N, as 


F=)S°\,n.@N, and J=|]X. (5.109) 


Since we intend to extend the matrix material to elastoplasticity in the next chapter, 
which relies on different mechanisms for the deviatoric and volumetric contributions, we 
introduce the isochoric components of the principal strains 


Aa = (DTPA = [Oy 7%, (5.110) 
b 


and define the isochoric deformation gradient as 


F=ŅX »,n,0 No. (5.111) 


Given that these deformation measurements adequately characterize a typical elastic ma- 
trix material, we now focus on the fiber material. 


Discussions on higher-order contributions for the constitutive modeling of composites with 
flexural resistance have previously been addressed, i.e. in Asmanoglo and Menzel |10, 11], 
using the framework as provided in the preliminary work of Spencer and Soldatos |114] 
and Soldatos |113]. Since the formulation is in terms of the reference configuration, the 
non-orthogonality of the fibers in the actual configuration has no effect as long as the 
reference configuration is orthogonal. 


A completely different but very straightforward and physically illustrative approach was 
presented by dell’Isola |31]. Therein a detailed investigation of a Piola homogenization of 
a discrete system of extensional and rotational 2D springs is investigated. We will derive 
the exact same deformation measures within a strain-gradient continuum framework and 
as a general extension to 3D, and show that it can be recast into the formulation of 
Asmanoglo and Menzel |10, 11]. 


As with the shells, we introduce two fiber threads oriented in the L and M directions, 
respectively, where L and M are constant orthogonal unit vector fields within the body 
in the reference configuration. We then derive the stretch of the respective fiber 


AL = || = |FL|| and Ay = [ml] = |F M|]. (5.112) 
With the normalized fiber directions 
~ l FL FM 
l= — and m= Z (5.113) 


à IFEI Au [FM] 
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we then define the change of the angle between both fiber threads as 


y =acos(l - m) — 2 


2 
= acos a ont) ae un 
= IFLIIFMI) 2 


Eventually the corresponding measure for the shear deformation of the fiber material 
yields 
$ = tan(y). (5.115) 


Within the context of generalized continua, specifically the strain-gradient theory, the 
fiber bending is described taking into account the gradients of the deformed fiber vectors, 
ie. VI = VFL and Vm = VFM. In particular, we consider 


VIL = \,ViIL+(VAL-L)i and VmM =\yVmM+(ViAm-M)m (5.116) 


which are projections of the fiber configuration gradients onto the initial fiber direction, cf. 
Asmanoglo and Menzel [10]. These expressions include terms related to stretch gradients 
of the fibers as well as fiber curvatures. We introduce the curvature measure for the fiber 
initially aligned in L-direction as 


i oe 
KL = E 

1 y 

= zV -l8 VAL)L (5.117) 
L 

1 FL FL 

= vrr- 9 (Tez 92) VF) 

[FLI? ( [PE] IFZ 


and for the fiber initially aligned in M-direction as 


1 
Km = —VmM 
ÀM 


1 z 
= yz (vm - m8 VAy)M (5.118) 
M 
1 FM FM 
= 5 VEM- o lo M) VF) M. 
FMI? ( [PMI NFM 


5.2.2 Variational formulation 
Energetic response 


The stored elastic energy density of the composite material is defined by the functional 


2 


tatia J) + gg (Ar, Anes du, Ku) (5.119) 
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where ¢ € [0,1] is the volume fraction of the matrix material. As a representative matrix 
material, we use the modified volumetric-isochoric decoupled Ogden material model 


Ü nat =y (FA, Ag, A3)) SE MS (JA, Ag, A3)) 


mat 


2 We ao ea) 


The parameters jz) and a» with b = {1,..., N} are related to the shear modulus and the 
parameters « and ß are related to the bulk modulus. Assuming that the fiber portion in 
both directions is identical, the corresponding contribution of the fiber material is defined 
by 
1 
Vey = 30 (AL - 1)? + Am — 17) +06” 


(5.121) 
+ 3 (KL "CKL + KM: C ky) 


where a and b are stiffness parameters related to stretch and shear of the fiber material. 
Moreover, the stiffness tensor related to fiber curvature is given as 


c=cyl@l+mem)+cnen with n=lxm. (5.122) 


Previous works such as Asmanoglo and Menzel |10], used a scalar bending parameter, 
resulting in an isodirectional bending stiffness and thus limited modeling capabilities. In 
the proposed tensorial argument c, the different stiffness parameters cy and cı can be 
interpreted as the in-plane and out-of-plane bending stiffness, respectively. These can 
in turn account for the geometric dependence, namely the respective plane moment of 
inertia of the fiber bundle, as detailled in Section 5.1.5. Note that we assume the same 
cross section for both fiber threads. 


Regarding the partial derivatives therein, we introduce first relations related to the Kirch- 
hoff stress T = Tmat + Tap as 


ray vol 
Tmat = T mat +7, T nat 


= N dev vol 
T mat ‚a + T mat a) Na ® Na 


Viat y Vina 
A E aye 


(5.123) 


and 


MT Oy, OF Du OF O6 OF Om, OF | Ory OF 


= Ç a OAL Men OAM Men 09 OV ar, OK, Vin > FT 
2 ’ 


(5.124) 
the higher-order stress of the fiber materials as 


Ba ==> = — (a i lal „Pak er) l (5.125) 


OK, OVE OKm OVE 
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Strong Form 


To identify and collect the internal contributions to the boundaries of the mechanical 
field, i.e. the resulting bending moments and normal stress contributions, we derive the 
internal virtual work as 


ow = A + ôV dV = f- : Vp + PB: VVdpdV. (5.126) 
Bo Bo 


The usage of the transformation 
7:V,69 = TF : Vip (5.127) 


along with integration by parts yields 


sy = [TAE GH Vip TFT): Vip +V: (Vö : P) dV. (5.128) 


Bo 
A second integration by parts related to the third term yields 


Bo 


In a last step, we apply the divergence theorem for the second and third term such that 
we obtain 


awe = [VO PTR N) apar [EFT -VP)N) 50 + BN : Vöpda. 


Bo OBo 
(5.130) 


Note that also contributions in tangential direction at the boundaries can be considered 
such that further integrations by parts incorporates the boundaries 0?Bp and 0*B, which 
represent curves and points, see e.g. Schulte et al. [106] and Javili et al. [67]. Assum- 
ing that the principle of virtual work 6W™ — 6W** = 0 is valid with respect to the 
corresponding functional spaces of admissible solution and test functions, the external 
contribution can be formulated as 


swe = | B.6p+ [T-pda+ | M:Vópaa, (5.131) 


T M 
Bo rg rà 


where B is a given body force per unit volume of the reference configuration. Moreover, 
T and M are a surface traction and a surface torque acting at the mechanical Neumann 
boundaries TT and TY. Eventually, we obtain the local form of the mechanical problem 
as 


V: (TF '-V-$)+B=0 (5.132) 
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1. Stress equilibrium 
V-(7F-7-V-$)+B=0 (5.136) 


2. Kirchhoff stress 


ow 
T = Tmat + Thb, Tmat = 2 Aa oe Na B Na; Tfib = zo IE FF FT (5.137) 


3. Dirichlet and Neumann conditions 


p=p(X,t) onTf, (TE -V-BN=T(X,t) on rf (6.138) 
VEN =V¢(X,t)N on TY”, PN = M(X,t) on pu 
4. Initial conditions SO DIE, (5.139) 
Table 5.3: Strong formulation of the problem. 
supplemented by boundary conditions 
p=p on Tf 
VeN=VOEN on TJ? (5.133) 


(TF '-V-P)N=T on IG 
PN=M on TY 


with prescribed fields p and VPN at the mechanical Dirichlet boundaries Tŷ and me: 


As usual for fourth-order boundary value problems, the entire boundary is decomposed 
twice, ie. To = T UTF with re OTF = Ø and To = TV UTM with Ty? NT = 4. 
For details related to the enforcement of the gradient condition given in (5.133), see e.g. 
Schuß et al. [108]. The problem in the strong form is summarized in 5.3. 


Weak Form 


Based on the derivations concerning the mechanical field within the previous section, the 
set of admissible test functions related to dp and its space are defined as 


V? = {op € H?(Bo) |p = OonT¥, VPN = OonTy*}, (5.134) 


included within the Sobolev functional space of square integrable functions and derivatives 
H* with k > 0. Then, the weak form of the multifield problem reads 


[ser + vv: 5p Bav - | õp: Taa- | Yop: Maa =o. (5.135) 
Bo rf ro 


Note that we neglect inertia terms within the mechanical balance equation, i.e. we consider 
only quasi static problems. 
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5.2.3 Numerical examples 


The following examples are dedicated to the verification of the higher-order bending 
contributions of the fiber material. In particular, we investigate the in-plane bending 
behavior using a benchmark test from the previous section with the shell formulation, 
as well as the out-of-plane bending behavior using a 4-point bending test. Additional 
examples demonstrating the capabilities of the proposed formulation for endless fiber- 
reinforced polymers are explored in the subsequent section. 


In-plane bending test 


=< 
A 2 


7 


Figure 5.32: In-plane bending test. Problem setting. The lines illustrate the fiber 
structure. 


We consider the same in-plane bending benchmark test as in Section 5.1.5, where the 
left edge of a Kirchhoff-Love shell is clamped while the right edge is subjected to an 
external in-plane torque, chosen to match a reference analytical solution. To verify the 
proposed formulation in terms of in-plane bending stiffness parameterization, we take the 
shell deformation result, extrude it to the corresponding 3D geometry and calculate the 
energy. The plate is of size L x W x H = 10 [mm] x 1 [mm] x 0.5 mm and is discretized by 
8 x 2x1 quadratic B-spline based elements, see Figure 5.32. Furthermore, we assume that 
the plate consists of a single fiber bundle with a cross section of A = HW = 0.5 mm? and 
a tensile stiffness of Eg, = 79000 N/ mm?. The area moments of inertia of the fiber bundle 
with respect to the e3-axis and the ez-axis are given by Ie, = HW?/12 = 0.0417 mm? 
and Ie, = WH?/12 = 0.0104 mm*, respectively. Using these quantities we calibrate the 
bending stiffness parameters as cy = Egple,/A = 6583.3333N and c, = Egpl.,/A = 
2212N. 


In Figure 5.33, the strain energy density is depicted for both the Kirchhoff-Love shell 
formulation as well as the proposed higher-order continuum formulation. Therein, we can 
observe the same homogeneous distributions which verifies the calibration of the in-plane 
bending stiffness parameter cy. Note that the parameter c, does not contribute to the 
simulation result, but will be investigated within the next example. 


4-point bending test 


Next, the out-of-plane bending behavior of the fiber material is investigated using a 4- 
point bending test. Therefore, we consider again a rectangular geometry of size L x W x 
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Figure 5.33: In-plane bending test. Strain energy distribution of the Kirchhoff-Love 
shell formulation (left) and higher-gradient continuum formulation (right). 
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Figure 5.34: In-plane bending test. Problem setting. The lines illustrate the fiber 
structure. 


H = 125 [mm] x 25 [mm] x 0.5 mm discretized by 50 x 10 x 2 quadratic B-spline based 
elements. The bidirectional composite material has a matrix volume ratio of ¢ = 0.53 and 
the fibers are aligned in the Y = 0° configuration, see Figure 5.34 for the details on the fiber 
orientation. The four point bending test as shown in Figure 5.35 leads to a pure out-of- 
plane bending deformation of the structure. In particular, we prevent the displacement in 
upward direction for the outer support points and prescribe a displacement in downward 
direction for the inner contact points. Additionally, the left support point is horizontally 
fixed, whereas we allow sliding for the other contact points. The material setting of 
the matrix material reads y = 1630.4N/mm? and a, = 2 for the deviatoric part and 
k = 6250 N/mm? and 8 = —2 for the volumetric part, which corresponds to a Young's 
modulus of Ema = 4500 N/mm? and a Poisson's ratio of y = 0.38. 
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Figure 5.35: 4-point bending test. Boundary conditions of the 4-point bending test. 
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Two different settings of the fiber material properties are applied assuming a single layer 
of fibers over thickness direction. Firstly, we set the tensile stiffness of the fibers to 
a = Ex, = 79000 N/mm? and the bending stiffness to cı = ON. Secondly, we set the 
tensile stiffness of the fibers to zero and adjust the bending stiffness as c, = Eq, H? /12= 
1645.83 N. 


The applied bending stiffness of the continuum fiber model correlates to the out-of-plane 
bending stiffness for a shell model with the same height. As shown in the previous exam- 
ple, the proposed strain-gradient continuum formulation matches the contributions of a 
gradient shell formulation, provided that the stiffness is chosen properly. Thus, if we re- 
solve the thickness of sufficiently flat geometry with in the continuum model to obtain the 
same deformation as expected for the shell theory, a coincident bending behavior of the 
structure should result. Figure 5.36 shows the load deflection result for the investigated 
material settings. As expected, both results match in a good agreement, i.e. the ten- 
sion/compression behavior of the continuum fiber model in this bending example can be 
described by the bending terms themselves. This is an important and well known result, 
as strain-gradient contributions emanate from a length-scale dependent microstructure 
and if this microstructure is already resolved by the first-order continuum framework, the 
second-order contributions must be removed. 
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Figure 5.36: 4-point bending test. Force-displacement curves for bending tests. 


6 A Hybrid Phase-Field Model for 
Fracture in Fiber-Reinforced 
Polymers 


Based on the concepts established in the preceding chapters, we now develop a comprehen- 
sive multifield framework for thermomechanical simulation of fiber-reinforced composites 
undergoing large deformations, and taking into account the microstructural size effects 
of the embeded fibers as well as inelastic, damage, and ultimate failure mechanisms. 
In particular, the polymeric matrix material is assumed to undergo porous-ductile frac- 
ture, while continuously embedded fibers undergo brittle fracture, as is characteristical 
e.g. for glass fiber-reinforced thermoplastics. A hybrid phase-field approach is developed 
and applied in conjunction with a modified Gurson-Tvergaard-Needelman GTN strain- 
gradient plasticity model that accounts for temperature-dependent growth of voids at the 
microscale. The mechanical response of the emerging woven fabric microstructure leads 
to additional higher-order terms representing homogenized bending contributions of the 
fibers. Considering the previous chapters, the underlying concepts and kinematic frame- 
works are only briefly introduced, with the focus instead on the multifield constitutive 
formulation for fiber-reinforced polymers. Eventually, a series of representative numerical 
examples, using a prototypical roving glass fiber-reinforced polyamide-6 (PA6) composite, 
is used to investigate the algorithmic performance and modeling capabilities, i.e. different 
kinds and sequences of failure within long fiber-reinforced polymers. 


6.1 Configuration and kinematics 


We consider a fiber-reinforced composite as a three-dimensional continuum body, along 
with the fields and kinematic framework introduced in Chapter 5 and Chapter 3. As- 
suming that the fiber deformation is congruent with the matrix deformation, we consider 
the joint field mapping (X,t) and the material deformation gradient F = Vy(X,t) 
with its determinant J = det(F) > 0. Moreover, we presume the absolute temperature 
6(X,t) as a further common field representing the thermal state of the matrix as well as 
the fiber material. 


Regarding the damage behavior of the matrix material, the variables introduced in Chap- 
ter 3 describe the porous plasticity and ductile fracture of the matrix material. Namely, 
we consider the equivalent plastic strain a( X,t) its dual, the dissipative resistance force 
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rr(X,t), the plastic deformation map F?(X,t), the void volume fraction f(X,t), and 
the crack phase-field s(X, t). 


Note that the biphasic character of semi-crystalline polymers, such as the matrix ma- 
terial polyamide-6 (PA6) considered in the numerical examples, is often modeled by a 
multi-mechanism (MM) constitutive model |20, 109] . Therein, the amorphous and crys- 
talline phases are each characterized by their own inelastic strain field and porosity vari- 
able. Due to the significant influence of the fibers, we opted for a more straightforward 
approach for modeling the matrix material’s mean mechanical effects, using a unified 
porous-elastoplastic model with averaged fields and material parameters, as in [73]. Fur- 
thermore, in agreement with the literature for the PA6 under study |68], porosity is 
assumed to be caused only by void growth. Consequently, void nucleation are neglected 
which simplifies the evolution laws of the isotropic hardening. 


In addition, we propose a hybrid phase-field formulation for the fracture modeling of 
the fibers. Assuming the woven fiber reinforcement structure as in chapter ?, with the 
constant and orthogonal unit vector fields in the reference configuration of the two fiber 
filaments, L and M, we introduce the order parameters 


sL(X,t) : Bo xT > R with SU [0, 1] and SL >0 (6.1) 


and 


su (X, t) : Bo xT>R with SM € [0, 1] and ŠM > 0 (6.2) 


describing the crack phase-field of the fiber aligned in L-direction and M-direction. As 
for the matrix material, the values s = 0 refers to the undamaged and s = 1 to the fully 
ruptured state of the fiber material. 


The variables introduced above characterize a multi-field setting for the formulation of 
temperature-dependent micro- and macromechanical damage in fiber-reinforced compos- 
ites based on seven independent fields 


Ir [p,0, a, r?, 5,51, 5m], (6.3) 


the finite deformation map «q, the absolute temperature field 6, the equivalent plastic 
strain field a, the dissipative plastic resistance force rP, the crack phase-field s of the ma- 
trix material, and the dual crack phase-field [51 , 5u] of the fiber material. The Lagrangian 
plastic deformation map F°? and the void volume fraction f will be condensed within the 
balance equations. 


The deformation measures for the porous-plastic fracture-degenerated matrix material 
are exactly the same as those introduced in Chapter 3. After a threefold multiplicative 
decomposition, namely an elasto-plastic decomposition, a volumetric-deviatoric decompo- 
sition and a fracture in-& sensitive decomposition over the eigenvalues, we finally obtain 
the elastic fracture insensitive part of the isochoric deformation gradient and the Jacobian 
determinant as 


- 3 . Oey af TIa>i 
F"=Y Xn N, and =) 2 


6.4 
TDS else 0 
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where g = ag((1—s)?—(1—s)”) —2(1—s)?+3(1—s)? is the adjustable degradation function 


via the modeling parameter ag. Moreover, A$ = (A°)9) are distortional stretches, defined 


via the isochoric, elastic parts of the principal stretches A% = 1105) 312, and lastly, 
the vectors Nna and N, represent the principal directions of the left and right stretch 


tensors. 


Regarding the fiber material, we use the deformation measures derived in Chapter 5. 
Along with the classical contributions for woven fiber structures considering the fiber 
directions, as in the stretches of the respective fiber threads 


AL =||FL|| and Au = ||FM|| (6.5) 


and the shear deformation 


din (aos Se: = z) | (6.6) 


we additionally consider the curvature measures that take into account the deformed fiber 
vectors VL = VFL and Vm = VFM 


1 FL FL 
KL = ——— VFL o (gge r) yF) 6.7 
LS TEL ( frz]? ri] en 
and 
1 FM FM 
Ku = ——— Mm o (SA eM): VF) M. (6.8) 
“= TEMP ( [FM] ® rm] 


Next, we incorporate the hybrid crack phase-field into the fiber kinematics. Assuming 
that the fiber material is brittle compared to the matrix material and that fiber rupture 
requires a local tensile state, we formulate fracture insensitive parts of the stretches as 


a yd, ee) if A 1 > Amen) if A 1 
4, = 4 Au) REIT nad Rees gh) E N (6.9) 
AL else Am else 


where we use the same type of adjustable degradation function as for the matrix Y. = 
dg, ((1—.)*—(1—s.)”) —2(1—s,)°+3(1—s,)*. Note that the e indicates the respective fiber 
direction L and M. For the shear deformation we propose the following degradation 


$ = gi (sL)gm(sm)o, (6.10) 


where in case of single fiber rupture the shear is completely degraded even if the remaining 
fiber is undamaged. Finally, the fracture insensitive measures of the fiber curvatures take 
the straighforward form 


KL = gL(sL)kL and Km = gm ($m) km. (6.11) 
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6.2 Variational formulation 


Next, we propose the constitutive framework for thermomechanical damage in fiber- 
reinforced composites. To be specific, we introduce constitutive energetic and dissipative 
response functions based on the above definitions and derive the required relations and 
evolution laws to formulate the multifield variational problem. 


6.2.1 Energetic response 


The stored thermoelastic energy density of the composite material is defined by the func- 
tional 


e me Je 1 eA/\ ï la fs 
Y = CULE, Je, 0) + SS uel e, Jans}, Bn, 0) 
Da a ls a ro t 
= 6 (Wo lE", F, 0) + Wu0)) + SS (Wine, A 8, Rr, Rar, 0) + Un (0)) 


(6.12) 
where ¢ € [0,1] is the volume fraction of the matrix material. The elastic contribution 
to the stored energy function of the matrix material is decomposed into volumetric and 
deviatoric parts 
mat mat 


ye = we (Far, dé, X,5),0) Lye (Fas, Ae; Ae, 5), 9) (6.13) 


As a representative non-linear constitutive law, a modified Ogden material model with 
the associated strain energy density function 


te=) 2, a (= -1) (6.14) 


and 
werd _ qa (An) H (J) 1) a OS b) (Coe 1) (6.15) 
Y 
is used for the numerical examples. The parameters and ay with b = {1,...,N} 


are related to the shear modulus and the parameters k and ß are related to the bulk 
modulus. Moreover, 04 is a reference temperature and the parameters e and y are related 
to the thermal expansion coefficient. Assuming that the fiber portion in both directions 
is identical, the corresponding elastic contribution of the fiber material is defined by 


1 ~ 7 pe 
fib = -a AL 1 2 | ÀM 1 7 | b 2 
5 x P+ (1?) +59 | dee 
+ 3 (RL + CK, + KM i cK) + av(6 = 00) (Os = 1) + (Am == 1) 


where v denotes the thermal expansion coefficient and a, b and c are stiffness parameters 
related to stretch, shear and curvature of the fiber material, as introduced in Chapter 5. 
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Next, the purely thermal contributions to the stored energy of the matrix and the fiber 


material are defined by 
8 0 
Viat = Cmat | 0 — 9% — 0 ln 7 (6.17) 
0 


0 
Wo, = 2c4r (0 — bo — 0 ln (=) (6.18) 


0 
respectively. Therein, Cma and Cfp are constant parameters representing specific heat 
capacities of the respective material. 


The evolution of the stored thermoelastic energy is given by 


£v = ¢ ( a Àe + an es Piat “aul 
C e 
a , (6.19) 
A OW as Vr. Min; (Vir + Van) 
— F F 0 
2 (3 te ep Fr Oe 


Regarding the partial derivatives therein, we introduce first relations related to the Kirch- 
hoff stress T = Tmat + Tab as 


— ray vol 
Tmat = Tmat +7, Tmat 


_ =. dev vol 
T mat, a op Tmat =) Na ® Na 


(6.20) 
> o ( Umar „ae 
>>. | EEE pee 
and 
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(6.21) 
the higher-order stress of the fiber material as 


_1-€ (OVE, ðk, OWE, Oña 
Po == = OVE * Iku OVE)’ vw) 


the driving force of the respective crack phase-field as 


en COW 1 — al 
— _ z= — = 2 
i ¢ Os ” Hi 2 os, , Hm 2 Osm i 3) 
and the specific entropy as 
N = Nmat + Nib 
= nat T That) = 1— C OV ii, + Vin) . (6.24) 


00 2 00 
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Moreover, we introduce a dissipation function 


Dint = Vpmat Tmat * d? + Vf, 


mat 


H5 + Veer, (HLS. + Husm), (6.25) 


to account for a transfer of dissipated energy due to plastification and fracture into the 
thermal field, where d? denotes the Eulerian plastic rate of deformation tensor. The 
above relations are derived in a thermodynamically consistent manner by assuming that 
the dissipated energy is completely transfered into the thermal field, i.e. by setting Vpmat = 
Vias = Via = 1. Note, however, that the plastic dissipation factor Vpmas is typically chosen 
in the range of 85% to 95% in the context of thermoplasticity, see e.g. [111, 130, 75]. In 
addition, based on experimental observations it is reasonable to set fracture dissipation 
factors to Vina < 1 and Vea < 1, see the discussion related to an energy transfer into the 
thermal field in |39, 107] and the references therein. 


To model the plastic and fracture mechanical response, we introduce an auxiliary func- 
tional as 

= x Faa 
Y = c (Tala, Va, 0) + Y (6, V5,0)) + Eh (sx, Vsr sm, Van) (6.26) 


mat 2 


As introduced in Chapter 3, the plastic contribution Dat describes the response of 


isotropic strain-gradient hardening related to the matrix material, with the equivalent 
plastic strain « its gradient Va and the plastic length scale lp 


3 i 2 
Tala. Va,0) = | yla) da + y(0)2| Val. (6.27) 
0 


The isotropic local strain hardening function y(a, 0) developed specifically for the semi- 
crystalline polyamide-6 in the framework of a multi-mechanism model |20, 109] is adapted 
for a unified approach (monomechanism) as in |73]. Furthermore the function is extended 
to thermoplasticity following |111, 102, 28]. In particular, we obtain the saturation-type 
function 


y(a, 8) = yo(0) + yi(0)explwpra] + ya(0) A — expl=wp2a]), (6.28) 


with the three temperature-dependent material parameters yọ > 0, yı > 0 and yo > 0 
defined as 


yo(0) = Yoldrer)(1 = wro(O = Oret)), 
yi (9) = Y1 (Ore) (1 >= (0 =, Pret); (6.29) 
ya(0) = YalOrer) 1 — cra (0 — Gree): 


Therein, the initial yield stress yo + yı determines the threshold of the effective elastic 
response, ya(0)(1—exp[—wp2]) describes an initial hardening stage. yı (0)exp[w,1a] allows 
for the simulation of rheo-hardening, in which large stretches of fibrils leads to an abrupt 
increase of stress. Moreover, wpı and wp2 are saturation parameters and wto, wt and wrz 
are thermal hardening/softening parameters. 
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Next, we formulate fracture contributions for the matrix as well as the fiber material. 
As in the previous chapters, we approximate a sharp crack surface T', by a regularized 
functional 


A 1 
neg i, False, Vs.) dV with False Vs.) = 5 (82 + Vse), (6:30) 
fe 
Bo 


based on a specific crack regularization profile Y. defined per unit volume of the reference 
configuration and the fracture length scale l which controls the regularization. Con- 
cerning ductile fracture of the matrix material, we require that lp > lẹ such that the 
regularized crack zone lies inside of the plastic zone. Using the regularization given in 
(6.30), the approximated fracture energy of the composite material reads 


A 1- > x 
wi I Jess Vs) + a (Ge, IL (SL; VSL) + Jem YM (SM; Vsm)) dV. (6.31) 


Here, g., denotes the Griffith-type critical energy density required to create fracture within 
the respective material. For the matrix material, the critical energy density is decomposed 
additively into elastic and plastic contributions as g.(@) = gep + Jee exp(—wra), with the 
modeling parameter wf. Summarized, the phase-field fracture contributions are given in 
terms of crack density functions as 


a 


U. == gela Yes, Vs) 


ae (6.32) 
= Sk + ve) 


and 2 
Vip = Ge, Yu (SL, VSL) + Goya (Sm, Vsm) 


_ Jer 7.2 2 2 Gent 7.2 2 2 
= Flot + la Iau?) + Sesh + R Vs?) 


Eventually, we obtain the dissipative resistance forces of the plastic field and the respective 
crack phase-field via the variational derivatives of W with respect to a and 5, as 


(6.33) 


rP = Côa DP ae = C(Oa Vr, — V (Ova Dat) (6.34) 


and 


rí = ô T =0, Y — V - (Op, V). (6.35) 


6.2.2 Dissipative response 


The composite’s dissipative response and evolution equations, which are very similar 
to the structure developed in Chapter 3, are now adapted and extended for fiber frac- 
ture. Regarding the porous elastoplastic behavior of the polyamide matrix material, 
multi-mechanism models |20, 109, 68] are commonly used in which local information such 
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as plastic strains, stresses, and damages are taken into account in each phase of semi- 
crystalline polymers via several GTN-type yield functions |57, 123, 96]. Since in the 
context of this work, we are only interested in the mean mechanical effects of the matrix 
material, a unified GTN-type yield function |73] is used, which implicitly defines averaged 
quantities, i.e. the effective scalar stress 4 := G(Omat, f) in terms of the Cauchy stress 
tensor Omat = Tmat/J and the void volume fraction f 


0? 
Ye (7, Omar, f) = Y + 2m fcosh KA —(1+(af)’) =0. (6.36) 


g2 


Here, oc = y 3/2 ||TAX/J|| denotes the von Mises equivalent stress, p = 4tr[Tmat/J] the 
pressure and q1/2 are fitting parameters. Note that for q, = 0 the influence of the pressure 
and the void volume fraction vanishes, i.e. © = 0. With the effective stress o and the 
dissipative resistance force rP we define the plastic yield function as 


PP (G(Omat, f) rP) =0 —rP. (6.37) 


As standard practice for polymide [68], we neglect other influences such as void nucleation 
or void softening due to shear, and obtain the evolution form of the void growth f = 
(1 — f)tr[d?]. Following [85], the current void volume fraction is given in terms of the 
plastic deformation as 

1— fo 


f=1-— 


A plastic Lagrange multiplier AP is introduced to enforce the Karush-Kuhn-Tucker con- 
ditions 


(6.38) 


P>0, OP <0, APO =0. (6.39) 


For the incorporation of the fracture mechanical behavior, we define crack threshold 
functions as 


PH. -= ri) =H.—r, (6.40) 


where the energetic driving forces He are bounded by crack resistance forces rf dual to 
the crack phase-field variables $e. Similar to plasticity, we introduce fracture Lagrange 
multipliers Af to enforce the Karush-Kuhn-Tucker conditions of the respective crack phase- 
field 

A>0, B <0, Aoi =0. (6.41) 


Based on the concept of maximum dissipation and the set € = [Omat rP, H — rf, Hy — 
ri Hm — rip AP, Af, AL, Afi], we define an extended dissipation potential and obtain a 
constrained optimization problem as where the Lagrange multipliers AP and àf control 
the non-smooth evolution of plasticity and fracture, respectively. Then, the associated 
plastic evolution equations follows as 


OP AP O®P 
d? = XP d à= -—— .42 
A Br and & IF or (6.42) 
and the evolution equation of the respective crack phase-field as 
det 
i (6.43) 
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A penalty regularization of the Lagrange multipliers can be utilized as follows* 


1 1 1 1 
AP = — (or ot el) APS 05 20) and Af = (04520) 

Np ne NE, Mu 
(6.44) 


where np, np and nr, are material parameters which characterize the viscosity of plastifi- 
cation and crack propagation. Note that in the sense of continuum setting as defined in 
(6.12) and (6.26), the rates obtained in (6.42) and (6.43) are weighted by the respective 
volume fraction ¢ and (1 — ¢)/2, respectively. 


6.2.3 Heat conduction 


Assuming that the absolute temperature 0( X,t) as a common field representing the ther- 
mal state of the matrix as well as the fiber material, we introduce a relation for the 
Piola-Kirchhoff heat flux vector as 


Q(F 0,5) = —K(F,s)V0 (6.45) 


which is known as Duhamel’s law of heat conduction. The thermal conductivity tensor is 
defined as 


K = (K(1 — 5) + K°"s) FF". (6.46) 
In case of fracture, the conduction degenerates locally such that we achieve a pure con- 
vection problem and the heat transfer depends mainly on the crack opening width of the 
matrix material. Here, we formulate the conductivity tensor K in terms of the phase-field 


parameter s. Moreover, K = (Kat + (1 — ¢)A gp is a average conductivity parameter 
related to the composite material and K“Y is a convection parameter. 


6.2.4 Coupled problem 


Based on the derivations concerning the mechanical field of the fiber-reinforced material 
just introduced, the set of admissible test functions related to Jl is given as 


SU = [59, 90, da, ôr”, ôs, 581, 554], (6.47) 


i.e. variations of the deformation, the absolute temperature, the equivalent plastic strain, 
the dual plastic resistance force, the crack phase-field of the matrix material and the 
variables of the dual crack phase-field of the fiber material, where their spaces are defined 


* The Macaulay brackets are defined by (x) = (x + |x])/2. 
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as 
V? = [Sp € H?(Bo) |öp = OonT%, VIPN = Oonly*}, 
V? = {60 € H!(B,)|60 = 0 onT¥}, 
V° = {da € H'(Bo)|da =0onIG}, 
v™ = {ör? € L’(Bo)}, (6.48) 
V’={ös € H'(B,)|6s =0 onT'}, 
Vv = !ös, € H'(Bo)| 05, = 0 onTr}, 
YM = {ősm € H'(Bo)| sm = 0 onm} 


included within the Sobolev functional space of square integrable functions and derivatives 
H* with k > 0. Then, the weak form of the coupled multifield problem reads 


[ee :r + VVP- ip: Bav - [sp Tas | Vip: Maa=0, 


Bo rz ry 
fwo — Dim — R) — V80 - Q dV — faas =0, 
Bo re 
J (Sa(Cy — r”) + Cyol, Vou - Va) dV =0, 
Bo 
.  Xp(&P)"? 

po (ma = EOE dV = 0, 
Bo 
[seme — ÓSXt (x — a) + xrCg-VOs-VsdV = 0, 
Bo f 

A 1-()ge 1- 
[rem — ÓSLX1, (m4 ( tl us.) IX Ook wer, - Vs, dV = 0, 
Bo ñ 

. 1 — ¿)ge 1- 
[rem — ÓSMX fu Ga _ as) + KG, Vis -VsydV =0, 

M 

Bo 


(6.49) 
where R is a given heat supply per unit volume of the reference configuration and Q is a 
heat supply across the thermal Neumann boundary re. For each other field, homogeneous 
Neumann boundary conditions are applied and appropriate thermal Dirichlet boundary 
conditions are formulated in terms of prescribed temperature 9, see Table 6.1. Note 
that we neglect inertia terms within the mechanical balance equation, i.e. we consider 
only quasi static problems. Additionally, internal conditions for the crack phase-field 
equations are given by 


s.=1 on Mch., (6.50) 
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ensuring that a fully broken state remains broken. The Karush-Kuhn Tucker conditions 
in (6.39) and (6.41), are evaluated by inserting local variables defined as 


1 for 9P>0 1 for $50 
Xp = and Xf =: A (6.51) 
0 otherwise 0 otherwise 


Using local variables xr, in (6.49), we demand $. > 0 for thermodynamical consistency, 
avoiding a transfer of dissipated energy back into the mechanical field. This prevents 
healing effects, which may be taken into account as well. We can also set xf, = 1 and 
restrict only the fully broken state, i.e. we allow for healing until the respective crack 
phase-field reaches the value one. 


6.3 Spatial discretization 


Due to the inclusion of curvature contributions in the fiber material, the mechanical part 
of the variational problem requires approximation functions, which are globally at least 
C*- continuous to satisfy yp", dp" € H?(BR). Therefore we apply the isogeometric analysis 
approach introduced in Section 2.5.2. Inserting the corresponding approximations and 
their variations ( yg, dp, sh, dsb, 0, 50%, at, da, ret, Sr?) into (6.49) yields the semi- 
discrete set of coupled equations 


da: / VV, RA +P: VVR4 dV — F&A! = 0, 


(0) 


50.4 Y (7° RAR? Op — RIDE, — VR4Q") dV — Q™*4 | =0, 


(0) 


da; | Ni (Gy — Nir?) dV + Klo; | =0, (6.66) 


> 


p | arii DR) 
or; | Maj — | xpN TL dV| =0, 
Bo 
SeA puts - / RH dV + KAP, B| =0. 

Bo 


Therein, rt, BP, nè and HA? are semi-discrete versions of the Kirchhoff stress tensor, the 
higher-order stress tensor, the local entropy and the phase-field driving forces obtained 
via the partial derivatives of the semi-discrete stored energy density 


ge = poe jo AA RRA), (6.67) 
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. Stress equilibrium 


a=0onTy, 


se = l on Ta, 


Initial conditions 


V-(TFT-V-B+B=0 6.52) 
. Kirchhoff stress 
ows. 1-COvÉ 
T = Tmat + Thb, Tmat = 2 AS Dre Na®Na, Tb 2 £ ap FT 6.53) 
. Higher-order stress ms 1-Cavs, ah 
— 2 OVE : 
. Energy balance 
O79 +V-Q-—Din —R =O 6.55) 
. Entropy 
0 (Weiss E Phat) L= Ço (WE, + Vin) 
N = Nmat + hb, Mmat ¢ 90 > Nfib 2 8 
. Dissipation . . . 
Dint = VpmarTmat : AP + Vena HS + Yen, (HL SL + Hmm) 6.56) 
. Piola-Kirchhoff heat flux 
Q=-KV0, K =(K(1—-s)+K°™s) FF" 6.57) 
. Plastic strain OP 1 z 
d? AP = 0, AP = (DP) P, Omat = Tmat/J 6.58) 
OT mat Np 
. Equivalent plastic strain AP aor 
A =0 6.59) 
. Plastic resistance force > 
rP = Côa VP at 6.60) 
. Crack phase-field equations 
. op! . 1 : 
eek O ALS of 6.61 
O(a — 7h) > iy Fe) 
. Crack phase-field driving forces awe 
He = — 6.62) 
OSe 
. Fracture resistance forces i sa 
Te = 05, Y 6.63) 
. Dirichlet and Neumann conditions 
p= p(X,t) on Tf, (FT V-BN =T(X,t) on TT 
Ve = Ve(X,t) on TY”, PN = M(X,t) on TY 
6 = 6(X,t) on TÍ, -Q:-N=Qo0n Te (6.64) 


Va:N=00nTy" 
Vse: N =00nT0 


p(X,0)=Po0, H(X, 0) =v, 0(X,0)=00, a(X,0)=0, r?(X,0) =0, s.(X,0)=0 (6.65) 


Table 6.1: Strong formulation of the coupled problem 
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cf. (6.12)-(6.24). DŁ, and Q" are semi-discrete definitions of the dissipation density and 
heat flux, cf. (6.25), (6.45) and (6.46). Moreover, the semi-discrete external contributions 


in (6.66); and (6.66), are formulated as 
EO / RAB AV + / RAT dA + / MVR? dA 
Bo LE Ta 


and 


get N RAR AV + J RAQ dA. 


Bo apn 


The coefficients of the matrices in (6.66)3 and (6.66)4 take the form 


Ke = Cyl J VN“.VN’dV and MẸ =n, f N° N’ dV, 
Bo Bo 


whereas the matrices in (6.66); are given by 


Ma? =n, / RAR? av, 


Bo 
Kern > fax (RARE + VRî^ VR?) av, 
E 
= 
| = sn | xs (RR? ai i VRA A VR?) dV, 
24, Pl 
E 
= = 21; Ga Xt (RARP + Ewe $ VR?) dV. 
M 
Bo 


(6.68) 


(6.69) 


(6.70) 


(6.71) 


Eventually, the semi-discrete functions Y*, ®? and g* denote the local hardening, the 


plastic yield and the critical fracture energy density, cf. (6.28) and (6.37). 


6.4 Temporal discretization 


Following the framework established in Section 3.4 ,the semi-discrete coupled problem 
(6.66) is discretized in time to obtain a set of non-linear algebraic equations to be solved 
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via a Newton-Raphson method 


ax | [rial +1 VR ay + PAU VES dV — Foe => 0, 


604 = er 2 wR Ros n+l — RED n+ 1- VRQ!) dpe QR E 0, 
Bo 


ar N’ (Cyh — NIP 4) dV + Kt ajnii| = 0, 


Mi Qj, mo Cant nun 2 (on™ _ 
ôr? | M, = - f xem Fed = fos) dV P= 0, 
Bo n+1 n+1 


fun 


Se Se 
958.4 MPAA ae ¿Bm - | ren; AT dV + KA ae — 0. 
Bo 
(6.72) 
Therein, a full-discrete definition of the internal dissipation using small values for the 
plastic viscosity parameter np, is given by 


h ] h a Uam+1 = ASAn+1 T An 
Dine n+l: = Vomat el = a) ret N I +v Vimat HER At 
h ASL,A,n+1 — SL,A,n h ASM,A,n+1 — 5M,A,n 
+ Uan Gar At T HM ny P At 


(6.73) 


Note that we apply a staggered scheme for the solution of the multi-field problem, i.e. the 
displacement field along with the plastic and hardening fields {qA n+1; Qin41; Tim +b FP ak 
the crack phase-fields s. an+ı and the temperature field 04 +1 are solved successively. 


For the time integration of the plastic evolution equations, the exponential integration 
scheme introduced in the context of the return-mapping algorithm in Section 3.4 is for- 
mulated in terms of the eigenvalues as 


OKS 


OO mat,an-+1 


Venti = Àg trEXP [-AtAR Mant] with nani = A (6.74) 


where Omatanti = (Tanzi + Tota n+1)/Jn+1- Note that in contrast to standard von 
Mises plasticity Ny +1 A Nir and [nny || # 1, i.e. the plastic correction has to be performed 
by the Lagrange multiplier AB, as well as the components Na n+1 which can be obtained 
by solving the non-linear relations 


OO 1 


OOmat,an-+1 


— nany =0 (6.75) 


n 


oP — M41 =9 and 


via an internal Newton-Raphson iteration. 
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6.5 Numerical examples 


In this section, we investigate the accuracy and performance of the proposed formulation 
for endless fiber-reinforced polymers. After verifying the higher-order contributions of the 
fiber material using two simple bending tests in the previous Chapter, we now focus on a 
series of tensile tests that demonstrate the ability of the proposed hybrid phase-field model 
to investigate different failure mechanisms for a roving glass composite material depending 
on the fiber configuration. This study is complemented by thermal investigations of the 
damage behavior of the model and its influence on the ultimate failure. Without loss of 
generality, we use an neo-Hook’s model for the matrix material in all examples, i.e., we 
put b = 1 and a; = 2 in (6.14). 


In particular, we consider a flat specimen of size L x W x H = 125 [mm] x 25 [mm] x 
2 [mm]. Figure 6.1 and 6.4 show the geometry in the reference configuration along with 
the applied boundary conditions and fiber configurations. Dirichlet boundary conditions 
are applied to the outer regions of length 20 [mm]. More specifically, one flap is fixed 
and the other flap is moved at a displacement rate of 0.5 [mm]/[s] within a quasi-static 
simulation setting neglecting inertial effects. The computational mesh consists of 2432 
square NURBS elements. The material setting, representing a prototypical roving glass 
fiber-reinforced polyamide-6 (PA6) composite, is summarized in table 6.2. We assume a 
square cross-section of fibers with Ag, = 0.0025 mm? and obtain a bending stiffness of 
CL = C4 = Egy Agiy /12 = 16.46 [N]. 


0 


A 


Figure 6.1: Tensile Test (unidirectional). Problem setting. The lines illustrate the 
fiber structure. 


6.5.1 Unidirectional fiber reinforcement 


We begin by examining the behavior of a unidirectionally reinforced composite with fiber 
orientations of Y = [0°, 10°, 20°, 30°, 40°, 65°, 90°], see Figure 6.1. The load deflection 
results for isothermal simulations at 0 = 293K are shown in Figure 6.2. It denotes 
the crack initialization and final fracture of the fiber material by O and o, respectively. 
Also, crack initialization and final fracture of the matrix material are indicated by © and 
x, respectively. Figure 6.3 shows the crack phase-field results of the fiber material and 
matrix material along with the plastic strain field results for the marked points. Note that 
the black mark indicates conditions without fiber fracture, since the specimen is already 
completely fractured. 
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Figure 6.2: Tensile Test (unidirectional). Load deflection results for unidirectional 
fiber reinforcements with different orientations. 


For a fiber orientation of Y = 0°, the fibers account for most of the load transfer due 
to the different Young’s modulus and fracture abruptly in the center of the specimen ©. 
Subsequently, the matrix material undergoes plastification and ductile fracture due to 
an abrupt load redistribution @. Notably, the increased strain rates result in a strong 
viscoplastic behavior within the matrix material, which is regulated by the viscous regu- 
larization parameter. 


In the unidirectional 10° fiber configuration, the fibers begin to crack near the clamping 
zones ®, which is further driven by the bending contribution to the crack driving force. 
This process is slowed down by the strain hardening behavior of the matrix material O. 
In the © state, the fibers are fully ruptured and the matrix material starts to fracture in 
the same area. 


At fiber orientation Y = 20°, the brittle fracture of the fibers starts again near the clamping 
zones ©. However, more plastification and thus solidification of the matrix material O 
occurs, so that the fiber and matrix material undergo final rupture nearly at the same 
deformation state ©. 


Applying a fiber orientation of Y = 30°, direct load transfer between both boundaries by 
the fibers is not possible since fibers which are clamped at the lower end do not reach the 
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Figure 6.3: Tensile Test (unidirectional). Results of the fiber crack phase-field (first 
row), the plastic strain field (second row) and crack phase-field of the matrix 
material (third row). The results are shown for the different deformation 
states and fiber configurations marked in Figure 6.2. 
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Figure 6.4: Tensile Test (bidirectional). Problem setting. The lines illustrate the 
fiber structure. 


upper clamping zone. Hence, the load has to be transferred towards the matrix material, 
resulting in higher plastification O and ductile fracture at the center of the specimen ©. 
Note that the fibers begin to fracture only in small areas near the clamping zones O. 


For fiber orientations of Y = [40°, 65°, 90°], the fiber material does not fracture due to 
small loads acting in fiber direction 0-0. Instead, the matrix material undergoes suitable 
plastification and subsequently ductile fracture leading to failure orthogonal to the fiber 
orientations. Concerning the Y = 90° fiber configuration, the fibers controls the necking 
which can be observed by comparing the deformation with the results obtained for pure 
matrix material O. This is also evident in the marginally increased stiffness of the load 
deflection results prior to crack initiation, where the results are almost equal up to a 
displacement of u = 40 [mm]. 


6.5.2 Bidirectional fiber reinforcement 


Next, we examine the same tensile test using a bidirectionally reinforced material with 
fiber orientations of Y = [0°, 10°, 20°, 30°, 45°] as shown in Figure 6.4. The load- defor- 
mation results for isothermal simulations at 6 = 293 [K] are shown in Figure 6.5. Again, 
the crack initialization and final fracture of the fiber material are represented by O and 
o, respectively. The final rupture of the matrix material is indicated by x. The crack 
phase-field results of the fiber and matrix material as well as the plastic strain results are 
shown in Figure 6.6. Note that only the phase-field results of the fiber oriented in the Y 
direction are plotted. 


For fiber orientations of Y = [0°, 10°, 20°], the bidirectionally reinforced material exhibits 
similar behavior to its corresponding unidirectionally reinforced counterparts. The ad- 
ditional orthogonal fiber merely provides higher necking resistance 0-0. In addition, 
the orthogonal fiber configuration restricts relative motion between the respective fibers, 
resulting in somewhat stiffer material behavior. This must be explored experimentally, 
which is outside the limits of the current work. 


For fiber orientations of Y = [30°, 45°], this effect becomes more pronounced, as can be 
seen in Figure 6.5. As previously discussed, the additional orthogonal oriented fibers 
counteract the necking behavior due to the Poisson effect of the matrix material, causing 
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Figure 6.5: Tensile Test (bidirectional). Load deflection results for bidirectional, 
orthotropic fiber reinforcements with different orientations. 


fractures within the matrix material near the clamping zones rather than in the center of 
the specimen 0-6. 


6.5.3 Thermal investigation 


Finally, we analyse the temperature dependence of the proposed model. For this purpose, 
we repeat the tension test with a unidirectional fiber reinforcement as shown in Figure 
6.1 using a fiber orientation of Y = 30°. 


Figure 6.7 depicts the load deflection result for isothermal simulations using temperatures 
of 0 = [253, 273, 293] [K]. The corresponding crack phase-field results of the fiber and 
matrix material as well as the plastic strain results are shown for the final deformation 
step in Figure 6.8. As previously observed, the matrix material experiences significant 
plastic deformations for 9 = 293 [K], followed by fiber fracture in small areas near the 
clamping zones, and finally the matrix material experiences ductile fracture in the center 
of the specimen. Reduced temperatures increase the yield stress of the matrix material, 
resulting in higher elastic energy and thus earlier, less ductile fracture of the matrix 
material. Note that in isothermal simulations with 9 = [253, 273] [K], the matrix material 
fails prior to the occurrence of any fiber cracks. 
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Figure 6.6: Tensile Test (bidirectional). Results of the fiber crack phase-field (first 
row) as well as the plastic strain field (second row) and crack phase-field 
of the matrix material (third row). The results are shown for the different 
deformation states and fiber configurations marked in Figure 6.5 
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Figure 6.7: Thermal investigations. Load deflection results for isothermal simulations 


at different temperatures of 0 = [253, 273, 293] [K] and a fiber orientation of 
Y = 30°. 
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Figure 6.8: Thermal investigations. Results of isothermal simulations at different 
temperatures of 0 = [253, 273, 293] [K] (each from left to right) and a fiber 
orientation of Y = 30°. Results are shown for the fiber crack phase-field 
(first block), the plastic strain field (second block) and crack phase-field of 


the matrix material (third block) at the last deformation state marked in 
Figure 6.7 


134 6 A Hybrid Phase-Field Model for Fracture in Fiber-Reinforced Polymers 


Elastic parameters 


Shear modulus u 1630 MPa 

Shear exponent a 2 

Bulk modulus x 6250 MPa 

Bulk parameter 3 -2 

Matrix volume ratio ¢ 0.53 

Tensile stiffness a 79000 MPa 

Shear stiffness b 0MPa 

Bending stiffness [c], c] [16.46, 16.46] N 
Plastic parameters 

Yield stress [Yo(Orer), Y1 (Orer), Y2(Orer)] (22, 56.8, 30] MPa 
Saturation exponent [wp1, wpa] [1, 115] 

Thermal softening parameter [wto, wt1, wt2] (0.4, 0.4, 0.4] K7! 
Viscoplastic parameter np 5000 MPa -s 
Viscoplastic exponent np 1 

Plastic length scale lp 3.1 mm 

Initial void fraction fo 0.01 

Gurson fitting parameter [q1, q2] [3, 0.8] 
Phase-field fracture parameters 

Brittle critical fracture energy [Jc e; Jer; Jem] [500, 500, 500] kJ /m? 
Ductile critical fracture energy gc,p 50 kJ/m? 
Saturation exponent wr 3 

Fracture viscosity [nr, nt, Nim] 1, 1, 1]-1077 MPa-s 
Fracture length scale [Ir, Ir, , liu] 3.1, 3.1, 3.1] mm 
Degradation parameter [ag, Ag, , Ag] 0.001, 0.001, 0.001] 
Thermal parameters 

Specific heat capacity [Cmat, Can] 1860, 2080] kJ/(m* - K) 
Thermal expansion coefficient [e, v] 106, 5] - 1076 K7! 
Thermal expansion parameter y 1 

Conductivity K 0.25 W/(m- K) 
Convection Kony 0W/(m-K) 
Reference temperature Oyof 293 K 

Fracture & plastic dissipation factor [Vomar Vimar» “fav! (0.9, 0.9, 0.9] 


Table 6.2: Material setting of the fiber-reinforced composite (PA 6/Roving glass). 


7 Summary and outlook 


7.1 Summary 


The non-linear framework presented in this work allows for a comprehensive investiga- 
tion of damage and fracture in fiber-reinforced polymers. The combination of a second- 
gradient theory, a novel hybrid phase-field model and a temperature dependent GTN- 
type plasticity model provides a numerical framework which is able to describe different 
failure mechanisms in detail. This approach allows for improvements in the design of 
such composite materials since we are able to predict fiber and matrix failure and their 
sequence dependent on the fiber orientation. Moreover, due to the fully-coupled, ther- 
momechanical microstructural approach we can optimize the fiber orientation for specific 
loads and thermal states. The verification and adjustment of the corresponding material 
parameters allows us to identify the physical effects of the higher-order terms in this pro- 
totypical second-gradient material. Several numerical tests conducted within this work 
have demonstrate the capability of the proposed framework to investigate such a complex 
behavior including the growth of microvoids, plasification and necking, crack initiation 
and propagation within the composite material and its components, respectively. 


7.2 Outlook 


Based on the present contribution the following research projects are currently in progress 
or in application: 


e A framework for simulating thermoplastics with short fiber-reinforcements, which 
are easier and cheaper to produce than continuous fiber composites, is currently 
under development 


e A variationally consistent multidimensional coupling approach for fiber-reinforced 
polymers where the Youngs moduli vary widely is under investigation 


e Contact mechanics for second- and third-gradient materials based on the concept 
of mortar-type methods are in progress 


e A research project on the simulation of the cracking process in bar-reinforced ultra- 
high performance concrete (UHPFRC) is currently in application 
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